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Abstract. In this review we focus on the determination of phase diagrams by 
computer simulation with particular attention to the fluid-solid and solid-solid 
equilibria. The methodology to compute the free energy of solid phases will be 
discussed. In particular, the Einstein crystal and Einstein molecule methodologies are 
described in a comprehensive way. It is shown that both methodologies yield the same 
free energies and that free energies of solid phases present noticeable finite size effects. 
In fact this is the case for hard spheres in the solid phase. Finite size corrections can 
be introduced, although in an approximate way, to correct for the dependence of the 
free energy on the size of the system. The computation of free energies of solid phases 
can be extended to molecular fluids. The procedure to compute free energies of solid 
phases of water (ices) will be described in detail. The free energies of ices Ih, II, III, IV, 
V, VI, VII, VIII, IX, XI and XII will be presented for the SPC/E and TIP4P models 
of water. Initial coexistence points leading to the determination of the phase diagram 
of water for these two models will be provided. Other methods to estimate the melting 
point of a solid, as the direct fluid-solid coexistence or simulations of the free surface of 
the solid will be discussed. It will be shown that the melting points of ice Ih for several 
water models, obtained from free energy calculations, direct coexistence simulations 
and free surface simulations, agree within their statistical uncertainty. Phase diagram 
calculations can indeed help to improve potential models of molecular fluids. For 
instance, for water, the potential model TIP4P/2005 can be regarded as an improved 
version of TIP4P. Here we will review some recent work on the phase diagram of the 
simplest ionic model, the restricted primitive model. Although originally devised to 
describe ionic liquids, the model is becoming quite popular to describe the behaviour 
of charged colloids. Besides the possibility of obtaining fluid-solid equilibria for simple 
protein models will be discussed. In these primitive models, the protein is described 
by a spherical potential with certain anisotropic bonding sites (patchy sites). 
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1. Introduction 

One of the first findings of computer simulation was the discovery of a fluid-solid 
transition for a system of hard spheres [H [2]. At that time the idea of a solid phase 
without the presence of attractive forces was not easily accepted. It took some time 
to accept it, and it was definitively proved after the work of Hoover and Ree [3], in 
which the location of the transition was determined beyond any doubt and, even more 
recently, when it was experimentally found for colloidal systems Certainly the 
study of phase transitions has always been a hot topic within the area of computer 
simulation. However, fluid-fluid phase transitions (liquid immiscibility, vapor-liquid) 
have received by far more attention than fluid-solid equilibria [5]. The appearance of 
the Gibbs ensemble [HI [7] in the late eighties provoked an explosion of papers dealing 
with vapor-liquid equilibria. The method has been applied to determine vapor-solid 
equilibria [8j, but not for studying fluid-solid equilibria. 

An interesting approach to the problem of the fluid-solid equilibrium was provided 
by Wilding and coworkers, who, in 1997, proposed the phase switch Monte Carlo method 
[9],[10]. This method was first applied to the study of the free energy difference between 
fee and hep close packed structures of hard spheres P [TT]. Three years later. Wilding 
and Bruce showed that the method could be applied to obtain fluid-solid equilibrium, 
and the fluid-solid equilibria of hard spheres was determined for different system sizes 
[T2| [13] . Quite recently Wilding and coworkers and Errington independently illustrated 
how the method could also be applied to Lennard- Jones (LJ) particles [m [15]. In our 
view, the phase switch Monte Carlo is closely related to the "Gibbs ensemble method" 
for fluid-solid equilibria, because, as in the Gibbs ensemble method, phase equilibria 
is computed without free energy calculations. In the phase switch methodology, trial 
moves are introduced within the Monte Carlo program, in which configurations obtained 
from simulations of the liquid are tested for the solid phase and vice- versa (phase switch). 
For a certain thermodynamic condition (p and T) the relative probability of the system 
being in the liquid or solid phase is evaluated and that allows one to estimate free energy 
differences. In the phase switch methodology the system jumps suddenly from the liquid 
to the solid in just one step. This method has been reviewed recently by Bruce and 
Wilding |10J, and has proved to be quite successful for hard spheres and LJ systems. 
It is likely that the methodology can also been applied to molecular systems although 
results have not been presented so far. It is not obvious whether the methodology can 
be used to determine solid-solid equilibria in complex systems. 

Another alternative route has emerged in recent years. Grochola fT6] proposed to 
establish a thermodynamic path connecting the liquid with the solid phase. Of course 
phase transitions should not occur along the path. If this is the case then it is possible 
to compute the free energy difference between the two phases. It is fair to say that 
Lovett [17] was the first to suggest such a path although it was Grochola who developed 
into a practical way. The system goes from the liquid to the solid not in one step (as 
in the phase switching methodology) but in a gradual way. Several variations have 
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been proposed so far and the method has been apphed successfully to LJ, electrolytes, 
aluminium IT8| [T9] and also to molecular systems, such as benzene [20]. At this 
stage it is not obvious whether it could also be applied to study solid-solid equilibria. 

Another approach to get free energies of solids is to use lattice dynamic methods 
[2T] . By diagonalizing the quadratic form of the Hamiltonian the system may be 
transformed into a collection of independent harmonic oscillators for which the free 
energy is easily obtained. This procedure allows one to estimate free energies at low 
temperatures and fails for discontinuous potentials and when anharmonic contributions 
become important (close to the melting point). The method has been used by Tanaka 
et al. [221 EH] to get the melting point of several water models. 

In this review we focus on the determination of phase equilibria between two phases, 
where at least one of them is a solid phase. Therefore, the goal is not just fluid- 
solid equilibria but also solid-solid equilibria. Free energy calculations allow in fact the 
determination of the global phase diagram of a system (fluid-solid and solid-sohd). In 
this methodology the free energy is determined for the two coexisting phases and the 
coexistence point is obtained by imposing the conditions of equal pressure, temperature 
and chemical potential. Usually the chemical potential of the liquid is obtained via 
thermodynamic integration. Different methods are used to determine the chemical 
potential of the solid. In their pioneering work. Hoover and Ree used the so called 
cell occupancy method [211 [3]. this method each molecule is restricted to its Wigner- 
Seitz cell, and the solid is expanded up to low densities [3]. One of the problems of this 
method is the appearance of a phase transition in the integration path (from the solid 
to the gas). The method was also applied to the LJ system [251 l26j . 

In the year 1984, Frenkel and Ladd proposed an alternative method, the Einstein 
crystal method [27J. In this method, that has become the standard method for 
determining free energies of solids, the change in free energy from the real crystal to 
an ideal Einstein crystal (in which there are not intermolecular interactions and where 
each molecule vibrates around its lattice point via an harmonic potential) is computed. 
Since the free energy of the reference ideal Einstein crystal is known analytically, it is 
possible to compute the free energy of the solid. If the equation of state (EOS) and 
free energies of both phases are known it is then possible to determine the conditions 
for the equilibrium between the two phases. Repeating the calculation at different 
thermodynamic conditions then it is possible to determine the phase diagram for a 
certain potential model. This route has often been used in the past for a number 
of simple models including hard ellipsoids [28], the Gay Berne model |29j, the hard 
Gaussian [30], hard dumbbells [SI [32l [33l [34] , hard sphero cylinders [35l [36l [37] , diatomic 
LJ models [381 ISHl SO], quadrupolar hard dumbells [H], hard flexible chains [121 113] . 
linear rigid chains [HI HS], chiral systems [16], quantum hard spheres [17], primitive 
models of water [18], electrolytes [l9l[50], benzene [5ll[52], propane [53] and idealised 
models of colloidal particles [SH EH EH EOl EH EH]- Some of the main findings of 
this research (up to 2000) have been reviewed by Monson and Kofke [5]. Forty years 
after the first determination of fluid-solid equilibria (for hard spheres [2]) the number 
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of models for which it has been determined is still small. The situation is even worse if 
one considers studies of phase diagram calculations (including both fluid-solid and solid- 
solid calculations) for models describing real molecules. Then the number of considered 
systems is quite small, comprising nitrogen [59], alkanes [60], fuUerenes [61], ionic salts 
[621 [63] and just in the last years carbon [M], silicon [65], sihca [66], [67], and hydrates 
[68]. 

As can be seen, water was missing and this is surprising taking into account its 
importance as solvent and as the medium where life occurs 169J. Although water has 
been studied in thousands of simulation studies since the pioneering works of Barker and 
Watts [To] and Rahman and Stillinger |71j , the study of its phase diagram by computer 
simulation has not received much attention. Interest has focused mainly on the possible 
existence of a liquid-liquid equilibria [72l [TH [ZH [751 ESI EZl [78], [79] . The interest in the 
solid phases of water has been rather limited although one observes a clear revival in 
the last decade [801 [ED [821 [83l [82 [851 [86l [871 [88l [891 [^ 

llOOl llOll 11021 11031 1104j . The only attempt previous to our work to determine the phase 
diagram of water was performed by Baez and Clancy in 1995 |105[ I106j . Estimates of 
the melting point of TIP4P were provided by Tanaka et al. [22] [23] . Vlot et al. |107j . 
and Haymet et al. [108[ 1109] . Motivated by this our group has undertaken the task of 
determining the phase diagram for a number of water models jllOl llllj . The study of 
water revealed that phase diagram calculations are indeed feasible for molecular systems 
and that they constitute a severe test for potential models. It is clear that the phase 
diagram contains information about the intermolecular interactions |112l 11131 1114j . 

The determination of a phase diagram is not, in principle, a difficult task. However, 
it is cumbersome, and somewhat tricky. In this work we will illustrate the details 
leading to the determination of the phase diagram of water. They can indeed be useful 
for those interested in water and its phase diagram. But the described methodology 
can be applied to other substances/models as well. We believe that by describing the 
calculations for water we are also describing how to do it for any other type of molecule. 
Problems where the determination of the fluid-solid equilibria by molecular simulation 
can indeed bring new light are among others, the design of model potentials for water 
and other molecules [69l [IISl Ill6l [20], the study of nucleation [ml [IIHl [Il9l Il20l [T2T] 
(where the equilibrium conditions should be known in advance), the study of the fluid- 
solid equilibria in colloidal systems, and also the very interesting problem of protein 
crystallisation. Our goal here is to describe all the details to encourage the reader to 
implement phase diagram calculations (including at least one solid phase) either to gain 
new insight on appealing problems or to improve currently available potential models. 

2. Basic definitions 

For a pure substance, two phases (labelled as I and II) are in equilibrium when their 
pressures, temperatures and chemical potentials are equal. A phase diagram is just 
a plot (for instance in the p, T plane) of the coexistence points between the different 
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phases of the system (gas, hquid or sohd). In this paper we shall focus on determining 
the phase equilibria for rigid molecules. Two ensembles are particularly useful to study 
phase equilibria: the canonical ensemble (NVT) and the isobaric isothermal ensemble 
(NpT). In the canonical ensemble the Helmholtz free energy A is given by the following 
expression [1221 126] : 

A = -kBTln{Q{N,V,T)) = -fc^Tln J exp[-/?[/(ri, r^, ^^)]rfl...t/ivj (1) 

where (3 = l/lksT), U is the intermolecular energy of the whole system, q is the 
molecular partition function and di stands for dviduji, where dvi = dxidyidzi. The 
location of molecule i is given by the Cartesian coordinates Xi,yi,Zi of the reference 
point and a normalised set of angles defining the orientation of the molecule (uti). 
By normalised we mean that / duJi = 1. For instance a reasonable choice for Ui is 
uji = fli/Vn where fli are the Euler angles |123j defining the orientation of the molecule 
and Vq = / dfli = Svr^. For a non-linear molecule the partition function can be written 
as [122] : 

q = qt'qrqvQe (2) 



q 



{2'nmkBTV''^ 



{2TikBTfl^Vn{hhhY'' 
s'h^ 



-p|- exp{-(3hiyj/2) 



In the previous equation translational and rotational degrees of freedom are treated 
classically (except for the symmetry number s' and for the factor /i), and vibrational 
and electronic degrees of freedom are described by the quantum partition function. 
qv = qt/V is the translational partition function (divided by the volume), and qr, qv and 
qe are the rotational, vibrational and electronic partition functions, respectively. The 
rotational, vibrational and electronic partition functions are dimensionless. We shall 
assume that the rotational, vibrational and electronic partition functions are identical 
in two coexistence phases. For this reason their precise value does not affect phase 
equilibria and we shall simply assume that their value is one (we do not pretend to 
determine absolute free energies but rather phase equilibria). The first factor g^/ has 
units of inverse volume or inverse cubic length. It is usually denoted as the inverse of 
the cubic de Broglie wave length |122j (i.e. 1/A'^). Therefore in this work g^/ is given 
by: 

" A3 " {h'^/{2TimkBT)fl^ 

In the NpT ensemble the Gibbs free energy G can be obtained as G = 
—kBT\n{Q{N,p,T)) where Q{N,p,T) is given by : 

q^ f3p f f 
Q{N,p,T) = J exp{-/3pV)dV J exp[-f3U{uJi,s^, ..,sn, ujn;'H.)]V^ dsiduJi.dsNdujNi4:) 

where Sj stands for the coordinates of the reference point of molecule i in simulation box 
units. The conversion from simulation box units Sj to Cartesian coordinates rj can be 
performed via the H matrix = Hsj (the volume of the system is just the determinant 
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of the H matrix). When performing Monte Carlo (MC) simulations of solid phases it is 
important that changes in the shape of the simulation box are allowed (i.e., changes in 
H) . This is usually denoted as anisotropic NpT simulations. They were first introduced 
within MD simulations by Parrinello and Rahman [1241 I125j . and extended to MC 
simulations by Yashonath and Rao |126j . In anisotropic NpT Monte Carlo the elements 
of the H matrix undergo random displacements, and that provokes a change both in 
the volume of the system and in the shape of the simulation box. Further details about 
the methodology can be found elsewhere |126[ I127j . The use of the anisotropic version 
of the NpT ensemble is absolutely required to simulate solid phases. It guarantees that 
the shape of the simulation box (and therefore that of the unit cell of the solid) is the 
equilibrium one. It also guarantees that the solid is under hydrostatic pressure and free 
of stress (the pressure tensor will then be diagonal with the three components being 
identical to the thermodynamic pressure). 

3. Fluid-solid equilibrium from NpT simulations? 

A possible way to determine the fluid-solid equilibria is by performing simulations at 
constant pressure and coohng the liquid until it freezes. However, it is very difficult to 
observe in computer simulations the formation of a crystal and this is especially true 
for molecular fluids [821 [83l I128j . The nucleation of the solid is an activated process 
and it may be difficult to observe within the time scale of the simulation. In fact even 
in real experiments super-cooled liquids are often found [128j . The other possibility 
is to heat the solid until it melts. Experimentally, when a solid is heated at constant 
pressure it always melts at the melting temperature (with only a few exceptions to this 
rule). In fact Bridgman [129[ stated in 1912: "It is impossible to superheat a crystalline 
phase with respect to the liquid". Unfortunately in computer simulations (in contrast 
to real experiments) one may superheat the solid before it melts. This is well known for 
hard spheres [130j (with pressure being the thermodynamic variable in question) and 
for Lennard- Jones particles [131j . The same is true for other systems such as water. 
For water models it has been found that in NpT runs ices melt at a temperature about 
90 K above the equilibrium melting point [13211133] . Similar results have been obtained 
for nitromethane [134j or NaCl In that respect NpT simulations provide an upper 
limit of the melting point. Introducing defects within the solid reduces considerably the 
amount of superheating [135[ I136[ . 

The difference between the results of NpT simulations and those found in 
experiments (summarised in the Bridgman's statement) is striking. The explanation to 
this puzzle is that in experiments melting occurs typically via heterogenous nucleation 
starting at interfaces (real solids do always have interfaces) whereas in NpT simulations 
it must occur via homogeneous nucleation (due to the absence of the interface), requiring 
a rather long time [1311 11371 1138j . Therefore new strategies must be proposed to obtain 
fluid-solid (or solid-solid) equilibria from simulations. The first possibility is to compute 
separately the free energy of the liquid and of the solid and determine the condition of 
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chemical equilibrium. The second is to introduce a liquid nucleus in contact with the 
solid (i.e., a seed) since that will eliminate the superheating. These two possibilities will 
be discussed in this paper. 

4. Thermodynamic integration: a general scheme to obtain free energies. 

In thermodynamic integration the free energy difference between two states/systems is 
obtained by integrating a certain thermodynamic function along the path connecting 
both states/systems [139] . The path connecting the two systems or states must be 
reversible. No first order phase transition should be found along the path. We shall 
distinguish (somewhat arbitrarily) two types of thermodynamic integration. In the first 
one the two states connected by the path possess the same Hamiltonian (i.e. interaction 
potential) and they just differ in the thermodynamic conditions (i.e. T, p...). We shall 
denote this type of integration as thermodynamic integration. In the second one, the 
thermodynamic conditions are the same for the final and initial conditions (i.e. the 
same p and T or same density and T), but the Hamiltonian (i.e, the intermolecular 
potential) will be different for the initial and final systems. This type of integration will 
be denoted as Hamiltonian thermodynamic integration. 

4-1. Thermodynamic integration 

Assuming that the free energy at a certain thermodynamic state is known, the free 
energy at another thermodynamic state is determined by establishing a reversible path 
connecting both states. For a closed system, with a fixed number of particles, two 
thermodynamic variables are needed to determine the state of the system (for instance 
p and T or V and T). In practice, it is convenient to keep one of the thermodynamic 
variables constant while performing the integration. 

4-1.1. Keeping T constant (integration along isotherms) Once the Helmholtz free 
energy at a certain reference density pi = N/Vi is known, the free energy at another 
density p2 = N/V2 (T being the same in both cases), can be obtained as : 



The integrand can be obtained in a simple way from NpT runs, isotropic for the fluid and 
anisotropic for solid phases |124[ 11251 1126j (so that the equilibrium density is obtained 
for different pressures). 

4.1.2. Keeping p constant (integration along isobars) In this integration the 
temperature of the system is modified while keeping constant the value of the pressure. 
In this way the Gibbs free energy G is obtained for any temperature along the isobar 
starting from an initial known value. The working expression is : 



A{p2,T) _ A{p,,T) 




■P2 pi^p) 

n ksTp^ 



dp 



(5) 



NknT NksT 



G{T2,p) _ G{T,,p) 
NkBT2 NksTi 



H{T) 



dT 



(6) 
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where H is the enthalpy. In practice several NpT simulations (anisotropic NpT for 
solids) are performed at different temperatures and the integrand is determined from 
the simulations. 

4- 1.3. Keeping the density constant (integration along isochores) In this case the 
density is constant and the temperature is modified. The working expression is : 
V^) ^ AiT^,V) _ rT. U(T) 
NkBT2 NkeTi M NkeT'^ ^ ' 

For fluids the integrand is easily obtained from NVT simulations. Although equation 
(jTj) is quite useful for fluid phases, it is not so useful for sohd phases. The reason is 
that for solids the density should be constant along the integration but the shape of the 
simulation box should be not (except for cubic solids). In fact the equilibrium shape 
of the unit cell (simulation box shape) changes when the temperature is modified at 
constant density. 



4-2. Hamiltonian integration 

In this type of integration the Hamiltonian of the system changes between the initial 
(A = 0) and the final state (A = 1). This can be accomplished by introducing a coupling 
parameter (A) into the interaction energy of the system. The interaction energy becomes 
then a function of this coupling parameter (f/(A)). The free energy of the system will 
be a function not only of the thermodynamic variables but also of A: 

exp[-pU{X)]dl...dN 



A{N,V,T, A) = -keTln 



q_ 



(8) 



By performing the derivative with respect to A in equation ([8]) one obtains: 

dA{N,V,T,X) /dU{X)\ 



dX \ dX I ^yrX 



(9) 

By integrating this differential equation one obtains: 

A(iV,V^,T,A = l) = A(iV,V^,T,A = 0)+ /'"VM^\ dX. (10) 

This equation gives the difference in free energy between two states with the same 
temperature and density but with different Hamiltonian (intermolecular potential). A 
similar equation can be obtained within the isobaric-isothermal ensemble. In this case 
the difference in Gibbs free energy between two systems with the same temperature and 
pressure and different Hamiltonian is given by: 

G(Ar,p,T,A = l)=G(iV,p,T,A = 0)+ f'^' /^^) dX. (11) 

Jx=o \ dX I ^^^^^^^ 



CONTENTS 



11 



5. The machinery in action. I. Obtaining the free energy of the Hquid 
phase. 

Here we shall briefly describe three possibilities to obtain the free energy of the liquid 
phase (there are other possibilities). The three routes considered are : thermodynamic 
integration, hamiltonian thermodynamic integration and the Widom test particle 
method. 

5.1. Thermodynamic integration 

When the density of a fluid tends to zero, the particles are far apart so that 
intermolecular interactions are irrelevant, and the system tends to an ideal gas. 
Therefore the free energy of the real fluid at a certain density p and temperature is 
given by : 



where the first three terms on the right hand side represent the ideal gas contribution to 
the free energy (a logarithmic correction to the Stirling's approximation was included) 
and the last term is the residual part (a residual property is defined as the difference 
between that of the system and that of an ideal gas at the same temperature and 
density). To derive equation (JT2l) the rotational, vibrational and electronic contributions 
to the partition function, equation (2), were set to one. The integrand in equation f|T2l) 
tends at low densities to the second virial coefficient. The first term on the right hand 
side of equation (fT2l) is just a reduced density and of course is dimensionless (although its 
numerical value depends on the value of A). To avoid phase transitions the integration 
along the isotherm should be performed at supercritical temperatures. Once the density 
of the liquid is achieved, one can integrate along an isochore to low temperatures. 

5.2. Free energy of liquids by Hamiltonian thermodynamic integration 

Let us label as A the system for which the free energy is known in the fluid phase, and B 
the system for which the free energy is unknown. By introducing a coupling parameter 
one can change from the Hamiltonian of A to the Hamiltonian of B : 



where A is a parameter ranging from (system A) to 1 (system B). According to equation 
(fTOj) . the free energy difference between A and B is given by: 



where < Ub~Ua >n,v,t,x can be obtained by performing NVT simulations for a certain 
value of A. The value of the integral is then obtained numerically. 




(12) 



U{\) = (1 - \)Ua + XUb. 



(13) 




(14) 
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5. 3. Widom test particle method 

The chemical potential can be obtained by the procedure proposed by Widom in 1963 
[140] which yields : 

This formula states that the residual value of the chemical potential fi^'^'^ is just the 
average of the Boltzmann factor of the interaction energy [Utest ) of a test particle. 
Although this formula is quite useful, its practical implementation may be problematic 
when the density of the system is high (so that inserting a particle becomes difficult). 
This is especially true for molecular systems and even more dramatic for systems with 
important orientational dependence in the pair potential. This is the case of water, for 
which it is quite difficult to obtain reliable chemical potentials by using the test particle 
method [T4T]. 

6. The machinery in action. II. Free energy of solids. 

The Einstein crystal method was proposed by Frenkel and Ladd in the year 1984 [27] 
and, since then, it has become the standard method to compute the free energy of solids 
[321 SH ^M, EHl [Hal ESI M, MM, [63] . in this method an ideal Einstein crystal is used 
as the reference system to compute the free energy of a solid. An ideal Einstein crystal is 
a solid (the word ideal pointing out the absence of intermolecular interactions) in which 
the particles (atoms or molecules) are bounded to their lattice positions and orientations 
by an harmonic potential and in which there are not interparticle interactions. The free 
energy of an ideal Einstein crystal can be computed analytically for atomic solids and 
numerically for molecular solids. 

For practical reasons, that will be clarified later, it is convenient to use an Einstein 
crystal where a certain reference point of the whole crystal is fixed. Two choices are 
possible: 

• Fixing the position of the center of mass. That was the original choice of Frenkel 
and Ladd [27]. We shall denote the reference system as an ideal Einstein crystal 
with fixed center of mass and the technique will be referred to as the Einstein 
crystal approach. 

• The second choice is to fix the position of just one of the molecules of the system, for 
instance molecule 1. In this second case we shall denote the reference system as an 
ideal Einstein molecule with fixed molecule 1. This methodology has been proposed 
quite recently [145] and will be denoted as the Einstein molecule approach. 

Since the determination of free energy for solids is rather involved and not many 
examples can be found in the literature describing the details, we shall describe both 
methodologies in detail. Obviously both approaches are quite similar and indeed provide 
identical values of the free energy of the solid phase. 
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(a) 



Solid 




AA3 


Solid with fixed CM 




AAi+ AA2 


Ideal Einstein crystal 


with fixed CM 




CM 




Ein-id 



(b) 



Solid 



-kTln(V/K] 



Solid with one fixed 
particle 



AA"i+ AA 2 



Ideal Einstein molecule 
with one fixed particle 



+)mn(y/K) 



Ideal Einstein molecule 

A— A Ein-mol-id 



Figure 1. Thermodynamic path used in (a) the Einstein crystal method (Frenkel 
and Ladd [27j and Poison eta/. [142J ). and in (b) the Einstein molecule approach 



6.1. The Einstein crystal method 

In the Einstein crystal approach the reference system is an ideal Einstein crystal with 
fixed center of mass. Let us describe briefly the notation that will be used in this section. 
The superscript CM indicates that the center of mass is flxed. The subscript specifles 
the interactions present in the system. In particular, the subscript Bin — id stands 
for ideal Einstein crystal (without intermolecular interactions), the subscript Bin — sol 
means that both the harmonic springs and the intermolecular interactions are present 
and the subscript sol indicates that only the intermolecular interactions are present 
(without the harmonic springs). 

The whole path from the reference ideal Einstein crystal with flxed center of mass 
to the crystal of interest can be described as : 

A-sol = ^Ein-id^ [i^Ein-sol ~ ^Ein-id) + i^'^'sof ~ ^Ein-sol)]^ (^soi " ^^t^ ) ■ ( 16) 

Here ^^in-jd is the free energy of the reference system (i.e. the ideal Einstein crystal 
with fixed center of mass). The first step is the computation of the free energy difference 
between the ideal Einstein crystal and the interacting Einstein crystal both with center 
of mass fixed (A^^fi^^.i - Ag;^„,J. In the second step {A'^f - A<^f^^^J, the springs of 
the interacting Einstein crystal are gradually turned off to obtain the crystal of interest 
(both systems with fixed center of mass). Finally the solid with fixed center of mass is 
transformed into a solid with no fixed center of mass {Asoi — A'^J}^). Equation (fT6l) can 
be written in a more simple way as : 

Asoi = A%f^_,, + [AAi + + AA3. (17) 
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By comparing equations (fT7|) and (fT6|) the meaning of the terms AAi, AA2, A A3 is 
clarified. Basically obtaining Asoi is a four step process, since you need to obtain 
^%ii-id (step 0), AAi (step 1), AA2 (step 2), and AAg (step 3). This integration 
path is schematically shown in figure [T] (a). 



6.1.1. Step 0. Obtaining the free energy of the ideal Einstein crystal with fixed center of 
mass: A'^f^_^^. As mentioned above, an ideal Einstein crystal is a solid in which the 
molecules are bounded to their lattice positions and orientations by harmonic springs. 
We will focus on rigid non-linear molecular solids. Although the translational field is 
always applied in the same form, the expression of the orientational field depends on 
the geometry of the considered molecule. We shall describe here the procedure for a 
molecule with point group C2V, as for instance water. The appropriate expression of the 
orientational field for other geometries will be given later on. 
The energy of the ideal Einstein crystal is given by: 

U Ein—id ~ U Ein—id,t ~\~ U Ein,or 
N 



u 



Ein—id,t 



E[ri - Ti 



(18) 
(19) 



U 



Ein,or 



2v 



N 
1=1 



N 

E 



Aij^aSin^ (?/'a,i) + J^E,b 



(20) 



In the preceding equation r, represents the instantaneous location of the reference 
point of molecule z, and is the equilibrium position of this reference point of molecule 
i in the crystal (i.e. will fiuctuate along the simulation run but not). A possible 
choice for the reference point (which defines the location of the molecule) is the molecular 
center of mass. In fact the rotational partition function of the molecule Qr is computed 
by using the principal moments of inertia (/i, I2 and J3) with respects to a body frame 
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with origin at the center of mass of the molecule. One could also use the center of mass 
of the molecule as the reference point to compute configurational properties. However, 
it should be pointed out that configurational properties do not depend on the choice of 
the reference point. For this reason, to compute configurational properties, there is a 
certain degree of freedom in choosing the reference point of the molecule. For free energy 
calculations it is very convenient (the reasons will be clarified later) if the reference point 
is chosen so that all elements of symmetry pass through it. This requirement is satisfied 
by the center of mass, but it may also be satisfied by other points. For instance, in 
the case of water, all elements of symmetry pass through the oxygen atom so that its 
choice as reference point is also quite convenient. Alternatively, one could argue that 
the oxygen would become the center of mass of the molecule if the hydrogen atoms 
would have zero mass. Since the phase diagram of a water model does not depend on 
the masses of the atoms forming the molecule, setting the masses of the hydrogen atoms 
to zero would not affect the phase equilibria. In this work we shall use the oxygen 
as the reference point of the water molecule and this choice will not affect the phase 
equilibria (take the reason you prefer either because configurational properties do not 
depend on the choice of the reference point, or because there is a certain combination 
of masses of the atoms of the molecule that render the center of mass on the oxygen 
atom ). The term UEin-id,t in equation (fT9l) is a harmonic field that tends to keep the 
particles at their lattice positions (rjo), while UEin,or forces the particles to have the 
right orientation. A^;, A^; ^ and Ke^ are the coupling parameters of the springs (not 
to be confused with the thermal de Broglie wave length A). Notice that A^;^^ and A.E,b 
have energy units whereas Ag has units of energy over a squared length. The angles ifja^i 

— * 

and Tph^i are defined in terms of two unit vectors, a and 6, that specify the orientation 
of the molecule, ipa^i is the angle formed by the unit vector a of molecule i in a given 
configuration (oj) and the unit vector (ajo) of that molecule in the reference lattice. The 
angle iph^i is defined analogously but with vector h. The definition of vectors a and h 
for a rigid triatomic molecule is shown in Figure [2l This form of the orientational field 
(equation ( |20l) ) was used by Vega and Monson [18] to get the free energy of a primitive 
model of water [ 146 ] . The vector a is calculated as the subtraction of the bond vectors 
a = {h — h)/ I h — h |, while b = {h + h)/ \ h + h |- The angles ijja^i and 'i/jb,i can 
be obtained simply from the scalar product of vectors Si and dio (both of them being 
unitary vectors), and 6j and bio (both of them being unitary vectors) respectively as: 

ipa,i = arccos {Si ■ Sio) 



so that ipa,i and ipb,i will adopt values between and vr. Notice that in the orientational 
field along the b direction (see equation (120|) ). the angle ipt,! is divided by a factor of vr. 
In this way, this term {il'b,i/'^Y ^iso takes values between (when b is parallel to bio) 
and 1 (when bi and bio form an angle of vr radians), and both orientational fields have 
the same strength (the sin^i^ipa) field changes from when -i/^a = or -^/^q = vr to 1 when 




(21) 



'a — 



vr/2). 
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The partition function of the ideal Einstein crystal in the canonical ensemble (after 
integrating over the rotational momenta) is given by: 



N „2 



QEin-id = J^^T^iqrQvqe)^ J exp dpi...dpN J exp [- pUEin-id] dl. . .dN , (22) 

L i=l 

where pi = {pxi,Pyi,Pzi) represents the momentum of molecule i and di = dridui, 
being rj the position vector of the reference point of molecule i and Ui its normalised 
angular coordinates. Consistent with our choice for the fluid phase, qr, qv, qe will be set 
arbitrarily to one (they will be omitted in what follows). Now a subtle issue appears. 
In the Einstein crystal approach each molecule (via the reference point) is attached to 
a lattice point. One can compute the free energy for a solid where each molecule is 
attached to one and only one lattice point. However one should not forget that there 
are A^! possible permutations. Therefore, the true free energy of the system is that 
obtained for a certain field where each molecule is attached to one lattice site multiplied 
by the number of possible permutations (i.e., A^!). For this reason the partition function 
is : 

QEin-id = -r^ exp -f3Y.7r~ dpi...dpN exp [-(3UEin-id]dl...dN, (23) 

iX J ZTfli done permutation 

where the integral over coordinates is now computed for just one permutation (and 
hence the label one permutation in the integral over coordinates). The expression one 
permutation in equation fl2^ reminds that each molecule is attached (via UEin-id) to 
one and only one lattice point. Let us now impose mathematically the condition of fixed 
center of mass of the reference points. For water we are fixing the center of mass of the 
oxygen atoms (the O will act as the reference point of the molecule) rather than fixing 
the center of mass of the whole system (including the hydrogens). It is simpler and more 
convenient for molecular fluids to fix the center of mass of the reference points, rather 
than fixing the center of mass of all the atoms of the system. In the configurational 
space, the restriction implies that: 

R'Cm(i"i, r2...rAr) — R^^j = 

N 

^/i,(r,-ri,) = (24) 

i=l 

where, if the mass assigned to all reference points (one per molecule) is the same, 
fii = 1/N, i = 1, ...N. In the previous equation Hcm is the center of mass of the reference 
points of the system (there is one reference point per molecule) in an instantaneous 
configuration and R^^^ is the center of mass of the reference points of the system 
when the molecules stand on the lattice positions of the Einstein crystal field. Rcm is a 
function of the coordinates of the particles of the system whereas R^-m is one parameter. 
Due to thermal vibration, in general Rca/ will be different from Rc-m- The constraint 
given by equation flM|) means that from all possible configurations of the particles of 
the system only those satisfying Rcm = Rcm '^i^ be allowed. 

A comment is in order here. The value of the molecular mass does not affect the 
phase equilibria (i.e., the molecular mass is irrelevant to determine phase transitions). 



CONTENTS 17 

For instance for a LJ system, the triple point does not depend on the particular value 
of the mass of the system (for Ar and Kr the melting point is different but not because 
of their mass but for the different value of the parameters of the LJ potential). For 
this reason, for phase diagram calculations it is quite convenient to assign the same 
mass to all molecules of the system (regardless of whether this is true or not in real 
experiments). For instance for NaCl, it is possible to assign the same mass to Na and 
CI, without affecting the phase equilibria of the model. In fact we have used this strategy 
to determine its melting point |147j . Therefore the simple choice /Xj = i = 1, ...N 
(i.e., assigning the same mass to all particles of the system or similarly to all reference 
points ) can be used to determine phase equilibria without affecting the results. We 
strongly recommend this choice. Of course dynamic properties depend on the mass, but 
not phase equilibria which is the main focus of this paper. 

As a consequence of the centre of mass constraint, the space of momenta is 
constrained to 

N 

Ep« = (25) 

i=l 

The partition function of an ideal Einstein crystal with fixed center of mass Q%fn-id 
be written as: 

QEin-id ~ Q Ein,tQ Ein,or (26) 

Then the free energy is simply obtained as: 

aCM _ aCM , . 
^Ein-id ~ ^Ein,t ^ ^Em,or 

= -keT In Q'^fl^-kBTlnQEin,or (27) 

The orientational term QEin,or will be computed by evaluating numerically the 
following integral : 



Q 



Ein,or 



1 f 1^ 

— r / exp {-PuEin,or) sm6d4>d6d-f 
57r^ J 



(2^ 



where (j), 9 and 7 stand for the Euler angles defining the orientation of the molecule and 
UEin,or IS the orieutatioual Einstein field for just one molecule (see equation (l20l) ). We 
have chosen to use the definition of Gray and Gubbins of the Euler angles |123] . In the 
particular case of a molecule with C2V symmetry (for instance water) it reads: 



Q 



Ein,or 



\ 2N "1 ^ 

'ipb 



J exp I -(3AE,a sin^ (ijja) - (3AE,b i-^j j sin 9d(t)d9d'-i 



(29) 



Notice that ipa and ipi, are functions of the Euler angles. The integral given by equation 
fl28|) or equation fl29l) can be evaluated numerically (for instance using a Monte Carlo 
numerical integration methodology). An approximate analytical expression [IS] has 
been provided for which is valid in the limit of large coupling constants (A^; A^; b). 

M 
in,t 

N 

6(^Pi)dpi...dpN 



The translational term Q%^t is given by the following expression: 



^ n2 



1 2mj 
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exp 



N 



N 



{Yi - rioj)dri...dr 



N- 



(30) 



Notice that to simplify the notation we have dropped the subindex "one permutation" 
in the integration over coordinates since it is sufficiently clear that this is indeed the 
case when each molecule is attached by harmonic springs to just one lattice points. This 
integral (equation fl5Ul) ) can be solved analytically [142J (see Appendix A for the details) 
and the result is: 

3(Af-l)/2 



Q 



CM 
Ein,t 



yCM 



Ml 



{N) 



3/2 



(31) 



The factor P^^^ accounts for the contribution of the integral over the space of momenta. 
Its value is not given explicitly, because we will see later that it cancels out with another 
similar term. 



6.1.2. Step 1. Free energy change between an interacting Einstein crystal and a non- 
interacting Einstein crystal (both with fixed center of mass): evaluating AAi. The free 
energy difference between two arbitrary systems 1 and 2 is given by: 

A,-A = -ksT In (32) 

Multiplying and dividing the numerator of the integrand by the factor exp(— /??7i), it is 
obtained that: 

A2-Ai = -ksT In (exp [-/3{U2 - f/i)])i (33) 

where (exp [—P{U2 — Ui)])^ is an average over the configurations visited by the system 
1. Taking U2 = UEin-id + Usoi and Ui = UEin-id (being Usoi the intermolecular potential 
of the solid), the previous expression can be written: 

A'^Efn-sol - A^L^d = 'k^T lu (cxp [- ^,01)]) E^n-M ^ 

Therefore, the free energy change can be computed simply as the ensemble average of 
the factor exp [—f3{Usoi)] along a simulation of the ideal Einstein crystal with fixed center 
of mass. This average is evaluated in a NVT MC simulation [HI HH [571 EH] • Note that 
this calculation must be done with the center of mass fixed. Often it is not possible to 
evaluate the free energy change as expressed in equation flMl) . because the exponential 
exp{—pUsoi) takes values larger than those that can be handled by a computer. This 
problem can be avoided if the expression is rewritten in such a way that the exponent 
does not take large values, for example, adding and subtracting from the energy of the 
solid Usoi the constant lattice energy Uiatuce- 

AAi = A%f^_^^, - A%f^_,, = Uiatuce - ksT In (exp [-/3(f/,,, - Ui^tuce)]) E^n-id ■ (35) 

One of the parameters that needs to be fixed when implementing the Einstein crystal 
method is the value of the spring constant (we will choose A^; = A^;^^ = Ae,^). A 
convenient choice for A^; is one that guarantees a small value (of about 0.02NkBT) 
for the second term on the right hand side of equation (l35ll . When this is the case 
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AAi is quite close to the lattice energy Uiatuce, defined as the intermolecular energy of 
the system when the molecules stand on the positions and orientations of the external 
Einstein field. 



6.1.3. Step 2. Free energy change between the solid and the interacting Einstein crystal 
(both with fixed center of mass) : evaluating /S.A2. The free energy change between 
the solid and the interacting Einstein crystal (both with fixed center of mass) will be 
computed by Hamiltonian thermodynamic integration. The harmonic springs are turned 
off gradually, and the total potential energy can be given by: 

U{\) = XUsol + (1 - X){UEin-^d + Usol). (36) 

The parameter A is defined between and 1, so that when A = one has the Einstein 
solid and when A = 1 one obtains the solid of interest. The free energy change along 
this path will be given by: 

AA2 = A{N,V,T,X = l)~A{N,V,T,X = 0)= f'~U^^) dX 

JX=0 \ dX /jvy,T,A 
l'X=l 

= — {UEin-id) N,V,T,\d'^- (37) 

It is a good idea to use the same values for Ae, AE,a, -^E,b- Then the spring constant 
along the integration are given by AA^, AA^; ^ and AA^; & (being all of them equal). It is 
convenient to perform a change in the independent variable from A to AA^; so that the 
integral of equation (1371) can be rewritten as: 

AA, = A-- - A-Z-.. = - r ^^""-"^"'"'"'^ ci(AA,). (38) 

Jo Ae 

Since the integrand changes several orders of magnitude it is convenient to perform a 
new change of variable [271 1139j AA^; to ln(AA£; + c) where c is a constant: 

Jln(c) I\E 

The integrand is now a smooth function of the variable ln(AA£; + c). A value [27] of 
c = exp(3.5) provides a good estimate of the integral (although the optimum value 
of c may depend on the particular considered problem). The integral of this smoother 
function can be accurately computed using, for example, the Gauss-Legendre quadrature 
formula [148] . It is usual to use between ten to twenty points to evaluate the integral. 
Fixing the position of the center of mass avoids the quasi- divergence of the integrand of 
equation ( l38l) when the coupling parameter A tends to zero. Without this constraint, the 
integrand would increase sharply in this limit (although it would remain finite), making 
the evaluation of the integral (equation (!38|) ) numerically involved, and making the 
accurate evaluation of the integrand at low values of the coupling parameter somewhat 
difficult. For this reason, it is numerically convenient to avoid the translation of the 
crystal as a whole for low values of A and this is achieved, either by fixing the center of 
mass, as in the Einstein crystal technique or by fixing the position of one molecule of 
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the system as in the Einstein molecule approach to be described below. In Appendix B, 
the procedure to implement the somewhat unpleasant condition of fixed center of mass 
within a Monte Carlo simulation is described. This is important since the calculations 
leading to A^i and Ay42 should be done with the center of mass fixed. 

6.1.4- Step 3. Free energy change between an unconstrained solid and the solid with 
fixed center of mass: evaluating AAs. As we have seen before (see equation (1321) ) the 
free energy change between two systems can be obtained as: 



AAs = - Af/ = -ksT\n^ = ksT In (40) 



where Qsoi is given (after integration over rotational momenta) by : 

^sol - JV!/^3(7V-1) j 



TV ^ r ^ 



1 2mi 



N 

6C^Pi)dpi...dpN 

1=1 

N 



/ exp [-(3Usoi{ri, uJi-.-r^, u^)] 5(^ fi,{ri - Yio))dviduji...dvNdujN- (41) 

i=i 

and Qsoi is given by an expression similar to that of but without the delta functions 
(and with h^^ in the denominator instead of h^^^~^^ ). Notice that the factor N\ cancels 
out when computing the free energy change (it appears both in Qsoi and Q'^^ )■ The 
integration over the space of momenta of the unconstrained solid is simply the integral 
of a product of Gaussian functions, whose solution is (when all molecules have the same 
mass) : 



(a) (^2' 



\ J 

The integral over the space of momenta of the solid with fixed center of mass is equal 
to the integral of momenta of the ideal Einstein crystal with fixed center of mass which 
was denoted as P*^*^. Substituting the partition functions in equation (140|) . we arrive 
to the following expression: 

/ pCM\ 

AA, = Asoi-AZf = kBT\n' 



P J 

^^^rj, /exp [-/3[/^o;(ri,u;i...rjv,^jv)](5(X:,=i(l/A^)(r» - Yio))dYiUJi...dr]^uji^ 
^ /exp [-/5f/soi(ri,u;i, ..vn^ujn)] dridui...drNduJN 

The energy of a system is not modified if the system is translated (while keeping the 
relative orientation of the molecules). The mathematical consequence of that is that 
?7soi(ri, uji, ..r^r, lun) can be rewritten as Usoii^i, i'25 ^2, ■■■'"^'n^n) where r- = rj — ri. Let 
us locate the center of mass of the lattice point at the origin of the coordinates system 
so that ( J2{'^/N)rio = R^m = 0)- Let us perform a change of variables from ri, r2, ...yn 
to v'2, ...HcM where Hqm is the position of the center of mass of the reference points. 
The Jacobian of this transformation is N. With these changes one obtains for the second 
term on the right hand side of equation P3|) : 



j^^rp / exp ( (cji , r ^ , a;2 , r 3 , . . . r ^ , c^Af ) ) (5 (Rca/ ) Ndui rfr ^ c^c^ac^rg . . dr'^duNdRcM ^^^-^ 
J exp{-pUsoi{uJi, r'2, LJ2, rg, ...r^, ujN))NdLJidr'2dLJ2dr'3..dr'j^dLJNdRcM 
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After integrating with respect to ujiy'2UJ2---'^'n^n one obtains: 




(45) 



since the Dirac Delta is normahsed to one. Now there is a quite subtle issue. The 
integral in the denominator of equation (HSil is just the volume available to the center 
of mass. What is the value of this volume? An interesting comment pointed out 
explicitly by Wilding [111 [IHl [HI [lO] is that the translation of a crystal as a whole 
under periodical boundary conditions generates N permutations between the particles. 
This is illustrated in figlS] for a two dimensional model. When counting the number of 
possible configurations we used the value A^! when going from equation (122] to equation 
f l23p . Therefore, we counted all possible permutations. Therefore the integral in the 
denominator of equation fH5l) is the volume available to the center of mass, within 
one given permutation. This value is simply V/N . Using V instead of (V/N) in the 
denominator of equation (H5l) is incorrect if the value A^! was used to count the number 
of permutations. In this case certain permutations would be counted twice, the first 
time in the factor A^! and the second via the translation of the whole crystal (in the 
volume V). Therefore : 



As can be seen the expression for AA3 is general and does not depend on the particular 
form of the intermolecular potential Usoi- Notice that correct results would also be 
obtained if one uses V in the denominator of equation (H5l) (so that the center of mass 
moves in the whole simulation box) but uses (A^ ^1)' when counting the number of 
permutations (i.e. count all permutations between particles except those obtained via 
the translation of the whole crystal through the periodical boundary conditions). In 
equation (!23|) one then would obtain a term (A^ — 1)!/A^! which provides an 1/A^ factor 
that could be joined with the ln(l/\^) term of equation (H5l) to give a contribution 
—kT\\i{y/N) which is identical to that given in equation fj46l) . Thus will have a 
term of the form — fcT ln(V^/A^) if A^! permutations were included in A^^g_^j^ (as done 
by Poison et al. |142j . and described here) or will have a term of the form —kT\n{y) if 
(A^ — 1)! permutations were included in A%^g_^^. Both choices are possible and provide 
the same total free energy. However when presenting results it is important to state 
clearly the choice not to confuse the reader. A sentence like that could be useful: 

• All permutations were included in the reference ideal Einstein crystal. That would 
indicate that a term A^! was used, and therefore contains a term of the form 
-kT\xi{V/N) 

• All permutations except those obtained by translation of the crystal under 
periodical boundary conditions were included in the reference ideal Einstein crystal. 
That would indicate that a term (A^ — 1)! was used, and therefore AAs contains a 
term of the form —kT\n{y) 

However when presenting results we recommend to join A^f^g_^^ and AAa into a unique 
term since the sum of both terms is unique and does not depend on the choice of the 



AA3 = Asoi - a; 



CM 
sol 



keT ln(P^*7P) - \n{V/N) 



(46) 
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number of permutations included in A^^^_^^. It is fair to say that Wilding [12] was the 
first pointing out explictly that N permutations were generated by translation under 
periodical boundary conditions. This has been taken into account implicitly by Poison 
et al. [132], since they used a term of the form —kT\n{y/N) for A A3. 



6.1.5. Final expression The final expression of the free energy of the solid is : 

Asoi = {A%'^Lrd + A^s) + AAi + AA2 = Ao + AAi + AA2 (47) 

where we have defined Aq as Aq = A%f^_^^ + AAa. Taking into account all the 
contributions to the free energy, the free energy of a molecular solid can be computed 
using the following expression: 
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- — In (exp [-l3{Usol - U lattice)]) i 
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A=o \ NkeT 



dX (48) 



Ny,T,X 



Nkj^T N — ^"'-^ I- /-v-a"' ^ i-u.i.Lice J i I Ein-id 

Notice that P^^^ does not appear in the final expression (so that its value is irrelevant 
for free energy calculations). The first two terms in equation (l48l) correspond to Aq. 
The last two terms on equation fHSl) are AAi and ^A2 respectively. The argument 
of the logarithm in the first term on the right hand side (embraced by brackets) is 
adimensional. In fact it has three factors, the first factor having dimensions of L~^-^, 
the second factor having dimensions of L^*^^~^) and the last factor having dimensions 
of L^. In any computer program a unit of length / is selected. It is quite convenient 
to set the thermal de Broglie wave length to A = /, and this choice should be used for 
the solid (in equation (HHj) ) and for the liquid (in equation (1T2|) ). Then the volume of 
the simulation box V (in equation ( l48ll ) should be given in units and the value of the 
translational spring A^ should be given in Energy/(/^) units. Notice that assigning an 
arbitrary value to A affects the absolute value of the free energies but it does not affect 
the coexistence properties. 

An important final comment is in order. Free energy calculations are usually 
performed in the NVT ensemble (with temperature and density fixed). It is quite 
important that the shape of the simulation box used in free energy calculations 
corresponds to that adopted by the system at equilibrium. It is not valid to impose 
(for instance from experiment) the shape of the simulation box since that will give free 
energies higher than the correct ones (the equilibrium shape minimises the free energy 
of the system for a certain density). Rather one should first perform NpT anisotropic 
Monte Carlo simulations |124[ I125[ 1126] , and determine the shape at equilibrium of the 
simulation box at a certain p,T and Hamiltonian (the density will be obtained as an 
average of the run) and then to perform free energy calculations in the NVT ensemble 
using the density and equilibrium shape of the simulation box obtained from the NpT 
runs. This remark is important for solids belonging to any crystalline class but cubic. 
A convenient choice for the vectors rjo, fli.o and hi^ that define the Einstein crystal 
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field (equations (fTSI) . (fT^ and (120]) ). is to use the equilibrium positions (to determine 
Fjo) and orientations (to determine a^^o and hi^) of the molecules of the system. Other 
choices are also possible (for instance fields driving the molecules into configurations 
slightly distorted from the equilibrium one). However the choice of the equilibrium 
configuration has the advantage that a lower value of the external field A^; is needed to 
obtain reliable results. Obviously the free energy of the solid should not depend on the 
particular choice of the vectors Fjo, ai,o and hi^ that define the Einstein crystal field. 



6.2. The Einstein molecule approach 

Quite recently Vega and Noya have proposed |145j a slightly different version of the 
Frenkel Ladd method. The method has been denoted as the Einstein molecule approach. 
The idea behind the Einstein molecule approach is to fix the position of one molecule 
of the system (say molecule 1) instead of fixing the center of mass. More precisely, by 
fixing the position of molecule 1, we mean that we fix the position of its reference point. 
The molecule can still rotate as far as its reference point remains fixed. Therefore we 
fix the position of molecule 1 (as given by the reference point) but we do not fix the 
orientation of molecule 1. Of course for a simple fluid (HS, LJ) there are no orientational 
degrees of freedom so that in the Einstein molecule approach, atom 1 is fixed. Fixing 
the position of one molecule avoids the quasi-divergence of the integrand of equation 
( |38|) when the coupling parameter A tends to zero. The computational implementation 
of the method as well as the derivation of the main equations is rather simple. 



6.2.1. The ideal Einstein molecule: definition and free energy The partition function 
in the canonical ensemble (after integrating over the space of momenta) is given by the 
following expression: 

Q = ^^^^sl / exp[-(3U{ri,uJi, ...,rN,uJN)]driduJi...drNduJN (49) 

We shall assign qr,qv,qe the arbitrary value of one. This expression can be written in 
a more convenient way by exploiting the fact that the potential energy of the system 
U depends only on the relative positions of the particles, but not on their absolute 
positions, i.e., it is invariant under translations of the whole system (while keeping the 
orientations of all the molecules in the translation). We will perform a change of variables 
from (ri, r2, Fat) to (ri,r2' = r2 — ri,...,riv' = vn — ri)- Under periodic boundary 
conditions and the minimum image convention, this change of variables leaves the limits 
of the integrals unchanged, because the maximum distance between two particles in any 
of the three directions of the space is always less than the length of the simulation box. 
Therefore: 

Q= ^,^3jv J dri J exp [-(3U{uJi,r'^,uj2, ■■■,r',^,uJN)]duJidr'2duj2...dr'^dujM 
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The value of the integral k is independent of the value of ri and, therefore, we can 
integrate over ri: 

Q = j^''- (51) 

The whole partition function (/t) can be computed by multiplying the integral 
corresponding to one permutation (/t') by the number of possible permutations, which, 
for a given fixed position of particle 1, is equal to {N — 1)!. Therefore, the partition 
function can be written as: 




Figure 3. Left : schematic representation of the Einstein molecule, in which particle 
1 is fixed and acts as the carrier of the lattice. The movement of all the remaining 
particles is given relative to the position of particle 1. Right: Permutations generated 
through periodical boundary conditions by the motion of particle 1. 

Let us now define the ideal Einstein molecule. The ideal Einstein molecule is an 
ideal system (without intermolecular interactions) where the reference point of one of 
the molecules (e.g., molecule 1) does not vibrate and acts as reference, while the rest 
of the molecules of the system (i.e., molecules 2,3,. .N) vibrate around their equilibrium 
configurations (see figure [3] for a schematic representation). The reference point of 
molecule 1 is called the carrier, because this point transports the lattice. Notice that 
in the Einstein molecule, molecule 1 can undergo orientational vibrations, as far as its 
reference point remains in a fixed position (obviously for a simple fluid there is no such a 
rotation, and the carrier is just the position of atom 1). The lattice (crystal) is uniquely 
determined by the position of the carrier. The Einstein molecule can move as a whole, 
and this motion is represented by the motion of the carrier, which is able to move and 
occupy any position in the simulation box. The expression of the energy of the ideal 
Einstein molecule is: 

U Ein—mol—id U Eins—mol—id,t ~\~ UEin,or 
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UEin-mol-id,t = X! ~ ^io)^ (53) 

1=2 

Notice that the main difference with equation (fT9|) is the absence of an harmonic term 
for the reference point of molecule 1. The orientational part of the potential is identical 
in the Einstein crystal and in the Einstein molecule approach. The partition function 
of the ideal Einstein molecule can be obtained by performing the integral k' for this 
particular case and substituting the value in equation (!52l) . The translational integral 
is particularly simple since is just a set of 3(A^ — 1) oscillators. The orientational 
contribution is obtained as in the Einstein crystal approach. Therefore, the Helmholtz 
free energy AEin-moi-id of the ideal Einstein molecule is given by: 



N N N \ V J 2 \ NJ \ 71 J N 

In the case that the carrier molecule (molecule 1) is fixed, the free energy will be equal 
to the free energy of the ideal Einstein molecule plus a term kT ln{V/A^) (where the 
term V comes from the constraint on the position of molecule 1, and the term comes 
from the constraint on the momentum). 



6.2.2. Integration path and computation of the free energy in each step. In the ideal 
Einstein molecule approach, the free energy of a given solid will be computed from 
integration to the ideal Einstein molecule. This integral is performed in several steps, 
that are summarised in the scheme shown in figure [H First the ideal Einstein molecule is 
transformed into an ideal Einstein molecule with one molecule fixed (what we mean by 
particle fixed is that its reference point remains fixed). Then the ideal Einstein molecule 
with one particle fixed is transformed into the real solid with one particle fixed. In the 
last step this fixed particle is allowed to move to obtain the real solid. As it is shown 
in the scheme, the factor fcTln(V^/A'^) that appears as a result of fixing one molecule 
in the ideal Einstein molecule cancels out with the free energy contribution of allowing 
molecule 1 to move to recover the real solid. As a result, the free energy of a solid can 
be computed simply by adding to the free energy of the ideal Einstein molecule, the free 
energy change between an ideal Einstein molecule with one fixed atom and the solid 
with one fixed atom (given by A A* + AAg): 

Asoi = AEin-moi-id + AA* + AA*^ = A* + AA* + AA*^ (55) 

where the asterisk in AA^ and AA2 serves to remind us that the integral should be 
performed while keeping the position of the reference point of molecule 1 fixed (and Aq 
is just AEin-moi-id) ■ The computation of the free energy change between the solid and 
the ideal Einstein molecule keeping one particle fixed is completely analogous to the 
computation of the free energy change between the solid and the ideal Einstein crystal 
keeping the center of mass of the system fixed. As in the Einstein crystal method, this 
free energy change will be calculated in two steps, represented by the terms AA\ and 
AA^. In particular, in the first step (AA*) we will compute the free energy change 
between the interacting Einstein molecule with one fixed particle and the ideal Einstein 
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molecule with one fixed particle. This free energy change is evaluated using the same 
procedure as in the Einstein crystal method with the difference that, instead of fixing 
the center of mass, the position of molecule 1 is kept fixed: 

AA* = Ulattice - ksT In (exp [-(3{Usol - UlatUce)]) Ein-mol-id ■ (56) 

which is formally identical to equation fl35|) except for the fact that averages should be 
computed over the ideal Einstein molecule system, rather than over the ideal Einstein 
crystal, and that molecule 1 will be fixed instead of the center of mass. 

In the second step, the free energy change between the interacting Einstein molecule 
with one fixed particle and the solid with one fixed particle is computed (AAg). This 
will be done by slowly switching off the springs of the interacting Einstein molecule : 

f/(A) = XUsol + (1 — ^){U Ein-mol-id + Usol) (57) 

The parameter A is defined between and 1, so that when A = one has the interacting 
Einstein molecule, and when A = 1 one obtains the solid of interest (both with the 
position of molecule 1 fixed). Usoi is the potential of the system under consideration. 
This equation is equivalent to equation (!36|) for the Einstein crystal. The free energy 
change in this first step will be calculated from the following expression: 

AAl = - ^^^"'-^'^^^'^•^■" ^(AA^). (58) 
Jo Ae 

which is identical to equation fl38l) except for the replacement of UEin-id by UEin-moi-id- 
The asterisk indicates that the reference point of molecule 1 is fixed in the integration. 
Notice that this integral does not diverge at low values of A, because the translations of 
the system as a whole are prevented by fixing the reference point of molecule 1. 

6. 3. Calculations for the hard sphere solid 

Let us now present some results for the free energy of a fee solid of hard spheres at a 
density p* = 1.04086. We shall compute the free energy using both, the Einstein crystal 
methodology [271 I142j described extensively in this paper and the Einstein molecule 
approach. Results are presented in Table [TJ The first point to be noted is that AAi and 
/S.A\ (and AA2 and AAg ) are similar but not identical (reflecting the fact that it is not 
exactly the same fixing the center of mass as fixing molecule 1). However, the sum of 
all terms contributing to Asoi gives the same value, so that the estimated free energy is 
the same (within statistical errors) with both methodologies. Obviously the free energy 
of a well-defined state should not depend on the procedure chosen to compute it. Since 
AAi and AA^ are quite similar, and the free energy of the system must be the same 
computed by both routes (fixing the center of mass or fixing molecule 1), then AA2 
and ^A2 must differ in about 3 In(A^) / (2A^) which is the analytical difference between 
Aq and Aq. This is indeed the case as it can be seen in Table [TJ The third aspect 
to be considered from the results of Table [H is that the total free energy presents a 
strong size dependence. Notice that this is not a problem of the methodology chosen 
to compute the free energy, but it is an intrinsic property of the HS solid (and likely 
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of other solids as well). In other words, the free energy of solids presents important 
finite size effects. This is further illustrated in figure H] where the free energy is plotted 
as a function of The estimated value of A/ [NksT) in the thermodynamic limit 

from our results is 4.9590(2), which is in good agreement with the estimates of Poison 
et al. [US] (4.9589), Chang and Sandler tl49j (4.9591), Almarza [EO] (4.9589) and de 
Miguel et al. (4.9586) ^151j (all obtained from free energy calculations although with 
different implementations). Therefore, the value of the free energy of hard spheres in 
the thermodynamic limit for the density p* = 1.04086 seems to be firmly established. 
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Figure 4. Left. Free energies of HS in the fee solid for p* = 1.04086 as a function 
of system size (filled circles). The open circles represent the free energies after 
including the Frenkel Ladd finite size corrections (i.e., adding (2/A^)lnA^ to the 
free energies of the solid). For all values of N the free energies obtained here (black 
circles) are in excellent agreement with those reported by Poison et al. |142] and de 
Miguel et al. [151]. Right. Coexistence pressure of the fluid-solid equilibria of hard 
spheres as a function of the system size as obtained from free energy calculations 
(filled circles) or from phase switching simulations as reported by Wilding [12J and 
Errington |l4j (open triangles). 



A consequence of the strong N dependence of the solid free energy is that the 
coexistence pressure p* also presents a strong N dependence as illustrated in figure HI 
It is of interest to estimate the properties at coexistence in the thermodynamic limit. 
We found [US], p* = p/{kBT/a^) = 11.54(4), p* = p^a^ = 1.0372, p* = pia^ = 0.9387, 
and /i* = p,/{kBT) = 16.04. The coexistence pressure is in agreement with estimates by 
Frenkel and Smit [131] (11.567), Wilding [12] (11.50(9)), Speedy [IS2] (11.55(11)) and 
Davidchack and Laird [153J (11.55). The Hoover and Ree estimate (11.70) seems now 
to be a little bit high. The chemical potential at coexistence obtained here is consistent 
with the value reported by Sweatman [154j (p* = 15.99 — 16.08) obtained using the self 
referential methodology to compute fluid-solid equilibria. 

Although finite size effects are present both in fluid and solid phases, they seem to 
be more pronounced in the solid (probably due to the coupling between the periodical 
boundary conditions and the geometry of the solid). In principle one is interested in 
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properties of the system in the thermodynamic hmit rather than for a finite size system. 
To estimate free energies in the thermodynamic hmit one should repeat the free energy 
calculations for several system sizes and extrapolate to the thermodynamic limit. This is 
quite involved and time consuming. For this reason it is of practical interest to introduce 
finite size corrections (FSC) that allow the estimation (although in an approximate way) 
of large systems free energies, by performing simulations of small systems (something 
similar to the g(r)=l approximation |130j used to correct for the introduction of the 
cut-off). Several recipes have been proposed recently [145] . Here we shall describe one 
of them, namely the Frenkel-Ladd FSC. 

Table 1. Free energy of the fee hard sphere solid at a density p* = 1.04086. 
The value of the different terms that contribute to the free energy in the Einstein 
molecule and in the Einstein crystal methods, are also shown. All free energies are 
given in NksT units. The thermal de Broglie wave length was set to A = u, the 
hard sphere diameter. 
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0.0018 


-3.6866 
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-3.6819 
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6.4- Finite size corrections: the Frenkel-Ladd approach. 

In the original paper of 1984, Frenkel and Ladd (FL) provided an expression for the 
free energy of the solid {2/N) InN higher than the correct free energy. That was first 
pointed out by Poison et al. |142j . In Appendix C the reasons for the appearance of 
the extra term {2/N) InN will be described. Thus the FL free energy Afjj^ /{NksT) is 
given by: 

Af^^/{NkBT) = A,,i/{NkBT) + {2/N) IniV (59) 

Notice that the term (2/A^)ln(A^) tends to zero in the thermodynamic limit, and 
therefore the FL expression is valid in this limit. However, for finite systems the FL 
expression gives a free energy higher than the true free energy of the system. For a typical 
system size N =350 the difference between both values is on the order of O.OSA^fcBT"- 
The interesting issue is that the FL free energies although incorrect (for a certain value 
of N) are relatively close to the value of the free energy in the thermodynamic limit. 
This is illustrated in figure H] for the HS system. For this reason, one may simply view the 
FL expression as containing an approximate prescription for the finite size corrections, 
providing free energies closer to the thermodynamic limit than the correct free energies 
of the system of finite size. Other approximate expressions for the finite size corrections 
(FSC) have been proposed recently [145] . 
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6. 5. The symmetry of the orientational field in Einstein crystal calculations 

For molecular fluids, the choice of the orientational external field used within the 
Einstein crystal (or Einstein molecule) simulations should be done with care since this 
can be a source of methodological errors. The position of the molecule is given by the 
position of the reference point. Things become simpler if the reference point is chosen 
in such a way that all elements of symmetry of the molecule contain this point. For 
instance, a convenient choice for H2O, NH3, benzene and N2 are the O, the N, the center 
of the hexagon, and the geometrical center of the N2, respectively. Now two orthogonal 
unit vectors a and b are attached to the reference point, and these two vectors are 
sufficient to define the orientation of the molecule (in fact two degrees of freedom are 
needed to locate the orientation of a unit vector a, and just one degree of freedom to 
locate the vector 6, which is perpendicular to a). Therefore the unit vectors a and b 
are a useful way of defining the orientation of the molecule (the Euler angles could be 
used as well but it is more convenient to use a and b). For convenience the vector b 
is chosen along the principal symmetry axis of the molecule (the C„ with the highest 
value of n). When the molecule i stands on its equilibrium position and orientation in 
the crystal then the vectors and bi adopt the values Ojo and bio respectively. Thus 
the subindex will refer to the orientation of the molecule in the equilibrium lattice 
position. Let us denote as ipa,i the angle between the vectors cij and aio and ipb,i the 
angle between the vectors bi and bio in an instantaneous configuration. Which expression 
should be used for the orientational field? The translational part will always be given as 
in equation ( fT9l) for the Einstein crystal approach or as in equation ( l53l) for the Einstein 
molecule approach. The orientational part will be the same for the Einstein crystal or 
Einstein molecule approaches. We have already given a convenient expression for the 
orientational field of a molecule with point group (as, for example, water). Let us 
now give a convenient expression for other symmetries. For a molecule with a point 
group of type a convenient expression for U Ein-id,or is : 



N 

UEin,or ^ ^ 
1=1 



^E,a sin' 



2 




(60) 



2 

For a molecule with a point group of type Dnh a convenient expression for UEin-id,or is 



UEin,or = ^S,aSin^ \ J + ^E,bSm'^ {'ipb,!) ■ (61) 

For a molecule with point group Oh, a convenient expression [57] for UEin~id,or is: 

N 

UEin,or = [^E,a sin^ {lpa,i,niin) + ^E,b siu^ {lpb,i,min) ■ (62) 
1=1 

where ipb,i,min stands for the minimum angle between bio and the six vectors connecting 
the reference point of the molecule with the six octahedral atoms/sites and an analogous 
definition for 4'a,i,min- For a linear molecule only one vector (i.e. vector bi) is needed 



CONTENTS 



30 



and the applied field should be of the form, for a D^ch '■ 

N 



UEin,or = \^E,bSm'^ 



(63) 



i=l 

For a Coo^v molecule a convenient choice is : 
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(64) 



When performing a MC run, it is convenient to introduce two different types of moves, 
translations and rotations. In a translation move the molecule moves as a whole and 
there is no change in the orientation of the molecule. Only the change in the translational 
energy with respect to the reference Einstein crystal (or molecule) needs to be computed. 
In a second type of move the molecule is rotated in a random direction and angle with 
respect to an axis passing through the reference point of the molecule. Since the reference 
point does not change the position under such a rotation, only the orientational energy 
with respect to the reference Einstein crystal (or molecule needs to be computed). 

The choice of an orientational field adapted to the symmetry of the molecule as 
the ones proposed here is highly recommended. When this is done the energy with 
the external field is invariant to any of the symmetry operations of the molecule. 
Thus, a standard MC or MD program will provide correct values of the orientational 
contribution to the free energy. One interesting questions is: is it possible to use an 
external orientational field that does not reflect the symmetry of the molecule? The 
answer is : in principle, yes, but you should write a special MC or MD code for that 
purpose. Special moves should be added where the symmetry operations of the molecule 
are implemented. For instance for water, one should incorporate the C2 operation that 
exchanges the positions of the two H atoms. Of course the energy of the molecule 
with the rest of the system is not affected by this operation. However, the energy of the 
molecule with the external orientational field may change when the external orientational 
field does not reflect the symmetry of the molecule (see the interesting paper by Schroer 
and Monson [51] illustrating this problem for benzene). 

For this reason, it is by far more convenient and simpler to use an orientational 
field that respects the symmetry of the molecule (examples for Cnv, Dnh, Oh and linear 
molecules have been given here). This subtle issue of the symmetry of the orientational 
field may have been an important source of errors in free energy calculations for 
molecular fluids. Let us just finish by saying that although we found convenient to 
have the vectors a and b orthogonal, other choices (as far as they are not colinear) are 
also valid and correct. 

6. 6. Einstein crystal calculations for disordered systems 

Let us now discuss briefly the case of disordered solids. When implementing the 
Einstein crystal harmonic springs are incorporated to fix the position (as given by the 
reference point) and the orientation of the molecules of the system to the equilibrium 
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configuration. In a disordered solid, there may be many possible "equilibrium 
configurations" differing in a significant way (not just differing in the labelling of the 
molecules). Let us just give three examples: 

Plastic crystals. Molecules with an almost spherical shape tend to form plastic 
crystals when freezing. In these plastic crystals the reference points of the molecules 
form a true lattice but the other atoms of the molecule are able to rotate (either with a 
free or with a hindered rotation) around the reference points. In principle the Einstein 
crystal methodology described in the previous section can also be applied to determine 
free energies for plastic crystals [28l[3Tl[32l[33l[36l[37l[57]. In addition to the translational 
field, an orientational field is included forcing the molecules to adopt an orientationally 
ordered solid for large values of the orientational field. Some issues should be taken into 
account when evaluating the free energy for a plastic crystal phase. At low values of the 
orientational field very long runs should be performed to guarantee that the molecules 
are able to rotate. Many values (20-30) of the coupling parameter A^; should be used 
to compute the integral of equations ( |38l) or (l58ll (the orientational contribution to the 
integrand increases quite rapidly for low values of the orientational field). Finally the 
absence of phase transitions along the integration path should be checked ( the external 
field should lead the system from an orientationally disordered solid at low values of 
coupling parameter to an orientationally ordered solid for large values of the coupling 
parameter without undergoing any phase transition). 

Water. In the case of solid water (say ice Ih) while the oxygens are ordered (i.e., 
they form a lattice) the hydrogens are disordered. However, Bernal and Fowler |155] and 
Pauhng |156] suggested that configurations satisfying the so called Bernal and Fowler 
rules have the same statistical weight and that configurations violating the Bernal Fowler 
rules can be neglected. The correct estimate of the experimental residual entropy of ice 
at K by using these two assumptions was a major achievement. Therefore the free 
energy of ice is approximated by : 

^ ^B-^ l^(^^con/jgMratio7i) 

-kBT\a{VL) - ksTlnlAcon figuration) (65) 

where A^onfiguration IS the free energy (obtained via the Einstein crystal methodology 
for a certain configuration satisfying the Bernal Fowler rules) and Q is the degeneracy. 
Pauling estimates —kBTln{Q) as —kBTln{{3/2)^). Therefore for ices one computes 
the free energy for a certain configuration and then adds the Pauling contribution to 
the free energy. This entropy can also be computed numerically (see for instance the 
work by Berg and Yang |157j ). Notice that when MC or MD runs are performed for a 
certain configuration of ice satisfying the Bernal Fowler rules, the system remains in this 
configuration along the run. This is because the time required by the system to jump 
from a configuration to another (both satisfying the Bernal Fowler rules) is beyond the 
typical time of a simulation run. equation (!65l) is useful not only for ice but for other 
disordered solids as well. In fact it can be applied successfully [1581 11591 1321 11601 138] 
to tangent dimers, formed by two tangent spheres, where Nagle |161] has estimated fl, 
and for fully fiexible hard sphere chains [12] (where Flory and Huggins |162^ I163j have 



CONTENTS 



32 



estimated Vt). 

Partially disordered phases. In certain cases the system possess disorder but, still 
certain configurations are more likely than others. Getting the free energy in such a 
situation is especially difficult. Firstly it is important to sample the configurational space 
properly to obtain equilibrium configurations of the system representative of the partial 
disorder. Secondly these configurations will differ in statistical weight, so that it does 
not seem a good idea to perform Einstein crystal calculations for just one configuration, 
since its statistical weight is unknown. In this case thermodynamic integration can be a 
more adequate route. An example of a partially disordered phase is the fee disordered 
structure of the RPM model (see discussion on this later on). 

7. The machinery in action. III. Obtaining coexistence lines: the Gibbs 
Duhem integration. 

Once the free energy of the liquid and the solid has been obtained for a reference state 
it is relatively straightforward to perform thermodynamic integration to obtain it for 
other thermodynamic states and locate a coexistence point between two phases (in 
case where it exists). The Gibbs-Duhem integration allows the determination of the 
coexistence lines once an initial coexistence point is known. 

7.1. Gibbs Duhem integration 

In the year 1993 Kofke realized that the Clapeyron equation can be integrated to 
determine coexistence lines |164[ 11651 I166j . The Clapeyron equation between two 
coexistence phases (labelled as I and II) can be written as: 
dp su - Si hjj - hi 

1^ = = V ^^^^ 

dl vii — vi I [Vii — vi) 

where we use lower case for thermodynamic properties per particle. Since the difference 
in enthalpy and volume between two phases can be determined easily (at a certain 
T and p) the equation can be integrated numerically. When implementing the Gibbs 
Duhem integration one obtains the coexistence pressure for the selected temperatures 
(the temperature acting as the independent variable). This is quite convenient when 
the coexistence line does not present a large slope in the p — T plane. When the slope 
of the coexistence line is large within a. p — T representation then it may be more 
convenient to integrate the Clapeyron equation in a different way: ^ = In this 

case the coexistence temperatures are determined for a set of selected pressures (the 
pressure acting as the independent variable). A fourth order Runge-Kutta algorithm 
is quite useful to integrate the differential equation. It is important to stress that 
anisotropic NpT simulations should be used for the solid phase within Gibbs Duhem 
calculations. Isotropic NpT simulations could be used for fluid phases and for solids of 
cubic symmetry. 
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7.2. Hamiltonian Gibbs Duhem integration 

When a coupling parameter A is introduced within the expression of the potential 
energy of the system, then a set of Generalised Clapeyron equations can be derived 
[3T| I167[ 1168] . For two phases at coexistence: 

gi{T,p, A) = gn{T,p, A). (67) 

If the system is perturbed slightly while preserving the coexistence it must hold that: 

vjdp - sjdT + j ^-^ = viidp - sndT + {^^^ ^.^^ (68) 

the last terms appearing in equation (l68l) are due to the presence of the new intensive 
thermodynamic variable A. If A is constant when performing the perturbation then one 
recovers the traditional Clapeyron equation. If the pressure remains constant when the 
perturbation is performed (so that T and A are changed) then one obtains : 

dT ^ T[{dgu/dX) - jdgj/dX)] 

dX hjj — hj 

It is simple to show (within the NpT ensemble ) that dg/dX is nothing but |^ = 
{ ^^dx^ )jv xx^ which can be determined within an NpT simulation. The final working 
expression of the Generalised Clapeyron equation (for perturbations of T and A while 
keeping p constant) is : 

dT ^ T{< dun{X)/dX >n,p,t,x - < dui{X)/dX >n,p,t,x) ^^^^ 

dX hii — hi 

This Generalised Clapeyron equation can be integrated numerically yielding the change 
in coexistence temperature (at a certain pressure) due to a perturbation of the 
Hamiltonian of the system (i.e., of the potential energy). A similar expression can 
be obtained for the case in which the system is perturbed at constant T (changing the 
pressure and A ). In this case one obtains : 

dp _ < duniX)/dX >n,p,t,x - < dui{X)/dX >n,p,t,x .^^x 

dX vji — vj 

The change in the coexistence pressure (at a certain temperature) due to a change in the 
Hamiltonian of the system (potential energy) is then obtained. Equations ( ITTi) and ( ITOl) 
will be denoted as Hamiltonian Gibbs Duhem integration. Hamiltonian Gibbs Duhem 
integration is a very powerful technique since it allows one to analyse the influence of 
the parameters of the potential on the coexistence properties. It also allows one to 
change the parameters of the potential to improve phase diagram predictions. These 
two possible applications will be illustrated later on for the case of water. 

In the particular case in which A is used as a coupling parameter taking the system 
from a certain potential to another (by changing A from zero to one): 

f/(A) = XUb + (1 - X)Ua. (72) 

Then the generalised Clapeyron equations can be written as: 

dT _ T(< Ub - UA >N,p,T,X - <ub-ua >n,p,t,x) 



dX hji — hi 



(73) 
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dp 



<ub-ua >n,p,t,x - <ub-ua >N,p,T,X 



(74) 



vii - Vi 



where is the internal energy per molecule when the interaction between particles is 
described hj Ub (with a similar definition for ua)- If a coexistence point is known for 
the system with potential A, then it is possible to determine the coexistence conditions 
for the system with potential B (it is just sufficient to integrate the previous equations 
changing A from zero to one). In this way the task of determining the phase diagram 
of system B (unknown) from the phase diagram of system A (known) is simplified 
considerably. 

8. Coexistence by interfaces 
8.1. Direct fluid-solid coexistence 

In 1978 Ladd and Woodcock devised a method to obtain fluid-solid equilibria, without 
free energy calculations, the direct coexistence method |169[ 11701 1171j . In this method, 
the fluid and the solid phases are introduced into the simulation box, and simulations are 
performed (NVE MD) to achieve equilibrium between the two coexistence phases. The 
coexistence conditions can then be obtained easily. Although the initial results for LJ 
and inverse twelve power were not very successful (probably due to the small size of the 
systems and to the short length of the runs), the method is becoming more popular in the 
last few years. In fact it has been applied to simple fluids |172[ 11731 11741 11751 1153[ I176j . 
metals [1771 [HH [1791 [HQ] , sihcon [M], ionic systems [1821 [183], hard dumbells [TM] . 
nitromethane [I35] and water [inHl IM HHHl [Ml [1871 [M IM Il90l IM]. Two 
simulation boxes, having an equilibrated sohd and liquid respectively, are joined along 
the z axis (the direction perpendicular to the plane of the interface). That could 
generate overlapping at the interface, and this overlapping should be relaxed/removed. 
The coexistence conditions (i.e., pressure, temperature) will be independent on the 
plane selected for the interface, but the dynamic behavior (and of course the interface 
properties) will be different for different planes |192l 11931 1190j . 

The direct coexistence method can be implemented either within Molecular 
Dynamics or Monte Carlo simulations. Both are equally valid, although if dynamical 
properties are of interest (for instance crystal growth rates) then MD is the only choice. 
The direct coexistence method was flrstly used in the NVE ensemble , but other 
ensembles as NVT, NpH, NpT, Np^T can be applied. Each ensemble will have its 
advantages and disadvantages, and the election of one ensemble or another depends 
on the information that one wants to obtain. Broadly speaking there are two kinds 
of ensembles, those at which it is possible to reach equilibrium having two coexistence 
phases at equilibrium and those for which it is not possible to have two phases at 
equilibrium. Obviously for the study of interface properties only those ensembles that 
lead to equilibrium should be used. 

A. NVE ensemble 
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This is the simplest approach. The idea behind the method is that the system 
will evolve to the equilibrium temperature and pressure by moving the interface (so 
that either the amount of solid or the amount of liquid increases). If the system 
is above the melting temperature, the solid will melt provoking (at constant E) a 
decrease of temperature. If the system is below the melting temperature, the fluid 
will freeze provoking (at constant E) an increase of temperature. In NVE runs the 
initial configuration should not be too far from equilibrium to guarantee some portion of 
liquid and solid in the final configuration. The equilibrium temperature and pressure are 
obtained in NVE simulations at the end of the run. The knowledge of the coexistence 
pressure only at the end of the run is a serious problem. In fact the lattice parameter 
used in the xy plane (which remains fixed along the run) may not correspond to the 
equilibrium value for the solid at the coexistence pressure. In other words, stress was 
introduced, and that may affect the free energy of the solid and, therefore, the melting 
point. This can be adjusted by trial an error |172[ 1153] . 

B. NVT ensemble 

In the NVT ensemble the system also evolves to equilibrium by changing the 
relative amount of liquid and solid phases, in this case to adjust the densities and 
pressure to their coexistence values. One important difference with the NVE ensemble 
is that, in the NVT, the heat released or absorbed by the crystallization or the melting 
is immediately accommodated by the thermostat and, therefore, it is expected that the 
system will attain the equilibrium faster than in the NVE simulations. Actually, heat 
transfer is usually the determining rate in the crystallization or melting process, and it 
has been seen that the presence of a thermostat in the simulations leads to crystallization 
rates much higher than those measured experimentally |194j . However, as in the NVE 
case, the solid is not able to relax in the xy plane and, therefore, the system might be 
under some stress. 

C. NpH ensemble 

A less common approach is to perform the simulations in the isobaric-isoenthalpic 
NpH ensemble |195j . with anisotropic scaling, i.e., the three edges of the simulation box 
change independently (see, for example, Ref. |187] ). In this ensemble, the system will 
also attain the equilibrium, in this case by evolving towards the coexistence temperature. 
One advantage with respect to the previous ensembles is that now the fluctuations of the 
volume will allow the solid to relax, removing the presence of stress. Moreover, as the 
volume of the box is allowed to change, the system can adapt more easily to changes in 
the relative ratio of the amount of solid and liquid phases, especially when the densities 
of the solid and liquid phases are very different. The problem in this case is that it 
is not strictly correct to use simulations under constant pressure in the presence of an 
interface, because, due to the contribution of the interface, the normal and tangential 
components of the stress tensor are not equal. However, if the system is chosen to be 
very large in the direction perpendicular to the interface, it is expected that the error 
introduced by the presence of the interface will be small. Another disadvantage of this 
ensemble is that, as in the NVE, the transfer of heat is not very efficient and, therefore. 
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long simulations are needed to obtain the equilibrium. 

D. NpT ensemble 

It is possible to tackle both the problem of having stress in the solid and of slow 
heat transfer by performing simulations in the anisotropic NpT ensemble where each 
of the edges of the box changes independently. In this the volume is able to 

fluctuate, the solid can relax to equilibrium and, as the temperature is fixed, the transfer 
of heat will occur very rapidly. However, as in the NpH ensemble, the use of the NpT 
ensemble in the presence of an interface is not strictly correct, although, as mentioned 
before, it is expected that the error is small for a system sufficiently large along the 
z axis. One important difference of this ensemble with the previous ones is that, as 
both the pressure and temperature are set, it is not possible to have the interface at 
equilibrium. The procedure to determine the coexistence properties is as follows. At 
a given pressure, different simulations are performed at a few temperatures. If the 
temperature is above the melting temperature the solid will melt (i.e. the total energy 
of the system will increase) and on the contrary, if the temperature is below the melting 
temperature the fluid will freeze (i.e. the total energy will decrease). In this way, it is 
possible to establish a lower and upper limit to the melting temperature. 

E. NpzT ensemble 

We have mentioned before that, due to the contribution of the interface, the stress 
tensor has different normal and transversal components and, therefore, it is not correct 
to perform simulations under constant pressure in the presence of an interface. The 
correct way would be to allow the size of the box to change only along the axis normal 
to the surface, i.e., the z axis. We will call such ensemble Np^T. The procedure to 
determine the coexistence properties is completely analogous to the procedure followed 
in the NpT ensemble. The only difference is that now a new starting configuration 
with the corresponding bulk lattice parameter at pressure Pz and temperature T must 
be generated for each simulation at different temperature and/or pressure in order to 
avoid having stress in the solid. This was not necessary in the NpT ensemble, as the 
fluctuations along the x and y planes allowed the sohd to relax to equilibrium. 

Two issues that deserve special attention when implementing the direct coexistence 
method are the system size and the length of the simulation. As with regards to the 
system size, a typical simulation box could have 10 molecular diameters in the x and 
y direction and about 30 in the z direction. Accordingly, studies by direct coexistence 
used typically 1000-3000 molecules, and these sizes provide results relatively close to the 
thermodynamic limit |173[ I178[ I135i 11891 1190j . Besides large system sizes, extremely 
long simulations are also needed (10 millions of time steps or more may be needed in 
many systems). Systems without a thermostat {NVE, NpH) may require even longer 
runs, since heat transfer along the interface may be quite slow. 
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8.2. Estimating melting points by studying the free surface 

It is now commonly accepted that melting starts at the surface and, already at 
temperatures lower than the bulk melting point, solids exhibit a liquid-like layer at the 
surface [HSl [Ml [HH [HH [2001 [M [2021 [2031 [ffl. Thus for most substances a liquid 
layer is presented in the surface even at temperatures below the melting point. When the 
thickness of this layer diverges at the melting point this is denoted as surface melting. 
When the thickness of the liquid layer remains finite (or even zero) at the melting point 
this is denoted as incomplete surface melting [205[ I206[ 1207] . The thickness of the quasi- 
liquid layer for a certain T depends on the considered material and on the exposed plane 
(as labelled by the Miller indexes). When the size of the liquid layer is sufficiently large 
either because the system has surface melting or incomplete surface melting (with an 
significant thickness of the liquid layer) then it is not possible to superheat a solid. 
This is the reason why experimentally solids usually melt at the melting point (at least 
one plane of the crystals of the powder presents a large quasi-liquid layer provoking 
the melting of the whole sample). For instance for ice it is not possible to superheat 
the solid, except for a few nanoseconds [208]. Superheating of solids (over macroscopic 
times) has been found experimentally only for monocrystals when the exposed planes 
have no liquid layer at all. Since for most of substances a quasi-liquid layer will be 
present at the melting point it is expected that for most of materials, the melting will 
occur at the melting point (when having a free surface). That provides a remarkably 
simple methodology to estimate melting points (at zero pressure). NVT simulations of 
the solid exhibiting a free surface are performed at different temperatures. A convenient 
geometry is to locate a slab of solid in the center of an orthorhombic simulation box. 
The pressure will be essentially zero since no vapor was introduced in the simulation 
box, and besides the vapor pressure is typically so small that the sublimation of a 
molecule from the solid will be a rare event within the typical simulation times. The 
lattice parameter of the solid in the direction perpendicular to the interface should 
correspond to the equilibrium values at zero pressure for the studied temperature. At 
temperatures below the melting point a stable thin liquid layer will be formed in the 
surface. At temperatures above the melting point the solid will melt. The simulations 
should be rather long to allow the system to melt completely. In the case of water, the 
ice took about 10-20ns to melt in the presence of a free surface. This technique has 
been applied successfully to estimate the melting point of the LJ model, water |209j and 
nitro-methane [2031 1210] . Notice that the technique will fail, for instance for NaCl, a 
substance having no liquid layer on the free surface [183] . 

9. Consistency checks 

Evaluation of free energy and coexistence points requires more effort than performing 
simple NpT runs. Besides, the possibility of introducing errors in the calculations is 
relatively high. For this reason it is a good idea to introduce several tests to guarantee 
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the accuracy of the calculations. 

9.1. Thermodynamic consistency 

In a number of cases it is possible to determine the free energy for two different 
thermodynamic states. For instance for a solid the free energy at two different densities 
along an isotherm can be determined by using Einstein crystal calculations. The free 
energy difference between these two states (as obtained from free energy calculations) 
should be identical to that obtained by thermodynamic integration (integrating the 
EOS along the isotherm). We indeed recommend the implementation of this test before 
performing any further calculation. Failing in the test indicates an error in the free 
energy calculations. However, passing the test it is not a definite proof of the correctness 
of the free energy calculations. They could still be wrong by a constant (being the 
constant identical for the two thermodynamic states) . It is clear that other tests should 
also be introduced to guarantee the correctness of the calculations. 

9.2. Consistency in the melting point obtained from different routes 

The melting point obtained from free energy calculations should be similar to that 
obtained from direct coexistence simulations (where the fluid and the solid phases 
coexist) within the simulation box. Notice that only fluid-solid equilibria can be studied 
by direct coexistence (it is not obvious how to implement solid-solid equilibria by direct 
coexistence). This is indeed a useful test. An incorrect prediction of the free energy 
of the solid phase (the typical source of errors) will provoke an incorrect prediction of 
the fluid-sohd equilibria as compared to the estimate obtained from direct coexistence 
techniques. Important differences (above 2-3%) in the melting point estimated from free 
energy calculations and from direct coexistence are not acceptable. The reader may be 
surprised by the fact that we stated that both melting points should be similar instead 
of stating that they should be identical. In the thermodynamic limit (for very large 
system sizes) they should indeed become identical. However for finite systems some 
technical aspects may provoke small differences. There are at least two reasons: 

(a) System size effects. The fluid-solid equilibria may have a strong N dependence. 
For this reason the coexistence pressure/temperature obtained from free energy 
calculations for a system of N molecules will not be identical to that obtained from direct 
coexistence simulations obtained with N* molecules. The case of HS was illustrated in 
a previous section. In general the stability of the sohd increases as the system becomes 
smaller. Since typically N* > N, then the sohd will appear shghtly more stable in the 
free energy calculations than in the direct coexistence results (assuming that no FSC 
corrections were introduced). This is what one may expect when the only difference 
between both type of calculations is the size of the system. 

(b) Cut-off' effects. This is important when the cut-off used in free energy 
calculations is different from that used in direct coexistence simulations. The difference 
in the melting point may just be due to the fact that we are simulating two sightly 
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different potentials. Even if the potentials were truncated at the same distance in 
both methodologies there may be differences. For instance, in free energy calculations, 
the potential may be truncated at a distance Tc, but long range corrections can be 
incorporated in the calculations |130j . In direct coexistence simulations, the potential 
may have been truncated at the same Tc but in this technique it is difficult to incorporate 
long range corrections. Therefore certain differences in the melting point between 
free energy simulations and direct coexistence calculations can be due to a different 
implementation of the potential. 

Notice that effects a) and b) may occur simultaneously. For instance in the case 
of water, different number of particles were used in free energy calculations and in 
direct coexistence simulations (effect a). But also the truncation of the potential (LJ 
contribution) was different in the free energy simulations (where long range corrections 
were included) and in the direct coexistence simulations (where they were not included) . 
In any case for water models we found differences in the melting point as estimated from 
both techniques of about 1-2% [mlfTSQ] . 

9.3. Some useful tests involving Gibbs Duhem integration 

Some useful tests that can be performed when using the Gibbs-Duhem technique are: 

• The coexistence lines should be identical (within statistical uncertainty) when the 
integration is performed forward (say by increasing the T) and when it is performed 
backward (say by decreasing the T). 

• If free energy calculations were used to locate an initial coexistence point between 
phases I-II, I-III and II-III, then the three coexistence lines obtained with Gibbs 
Duhem integration should cross at a point (the triple point). 

• If the melting point of two models has been determined by free energy calculations 
then it is also possible to use Hamiltonian Gibbs Duhem calculations to check that 
the melting point of model B is obtained starting from the melting point of model 
A and viceversa. 

• If the melting point of models A and B is known, then Hamiltonian Gibbs Duhem 
integration could be performed to estimate the melting point of model C. Both 
integrations (one starting from A and the other starting from B) should provide 
the same estimate of the melting point of C. 

• If two coexistence points between phases I and II are known (either by free energy 
calculations or by direct coexistence) a Gibbs-Duhem integration starting from one 
of them should pass through the other one. 

9.4. Consistency checks at K 

At zero temperature the condition of chemical equilibrium (i.e., the equilibrium pressure 
Peg) between two phases, labelled as phase I and II, respectively, is given by: 

UliPeg, T = 0) + PegVliPeg, T = 0) = Un{Peg, T = 0) + PegVll{peg, T = 0) (75) 



CONTENTS 



40 



Hence, phase transitions between solid phases at zero temperature occur with zero 
enthalpy change. This is really useful since it means that phase transitions at K can 
be estimated without free energy calculations (just computing mechanical properties 
as densities and internal energies). By performing several NpT simulations where the 
temperature is reduced up to zero it is possible to obtain the EOS (and internal energy) 
of each solid phase at K. Then by equating the enthalpy, it is possible to locate phase 
transitions (at K). This can be used as a consistency check. By performing free energy 
calculations it is possible to locate the coexistence pressure between two phases (I and 
II) at a finite non zero temperature. Then, by performing Gibbs Duhem simulations 
it is possible to determine the coexistence line up to K. The coexistence pressure at 
K obtained from this long route (free energy calculations+Gibbs Duhem integration) 
should be identical to that obtained from the short route (estimating the properties of 
the system at K). This is again a severe consistency check. Although runs at K 
enable a check to be made on the consistency of phase diagram calculations, they do 
not allow by themselves to draw the phase diagram of a certain model. Gibbs Duhem 
simulations can not be initiated from a known coexistence point at K since both Ai7 
and T are null so that its ratio AS*, which within classical statistical thermodynamics 
is finite even at K, can not be determined. 

10. A worked example. The phase diagram of water for the TIP4P and 
SPC/E models. 

We shall now illustrate how the previously described methodology can be applied to 
determine the phase diagram for two popular water models: SPC/E [212] and TIP4P 
[213] . We believe that they illustrate quite well the typical difficulties found when 
determining by computer simulation free energies of solid phases and phase diagram 
calculations. The SPC/E and TIP4P models are presented in Table [2] (along with 
two other recently proposed models TIP4P/Ice [214] and TIP4P/2005 [215] ). In these 
models a LJ center is located on the O atom, and positive charges are located on the H 
atoms. The negative charge is located at a distance doM from the O along the H-O-H 
bisector. Let us now describe briefly some of the simulation details. 

Table 2. Potential parameters for several water potentials. Notice that the OH 
bond length and the HOH angle are different for the SPC/E and TIP4P models. 



Model 


e/k{K) 


a(A) 


te(e) 


C?OAf(A) 


SPC/E pT2] 


78.20 


3.1656 


0.4238 





TIP4P 


78.0 


3.154 


0.520 


0.150 


TIP4P/2UU5 [215] 


93.2 


3.1589 


0.5564 


0.1546 


TIP4P/Ice [HI] 


106.1 


3.1668 


0.5897 


0.1577 
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10.1. Simulation details 

In our Monte Carlo simulations, the LJ potential was truncated for all phases at 8.5 
A. Standard long range corrections to the LJ energy were added. The importance of 
an adequate treatment of the long range coulombic forces when dealing with water 
simulations has been pointed out in recent studies |216l 12171 12181 1219] and it is 
likely that this is even more crucial when considering solid phases. In this work, the 
Ewald summation technique has been employed for the calculation of the long range 
electrostatic forces. The real space contribution of the Ewald sum was also truncated at 
8.5 A. The screening parameter and the number of vectors of reciprocal space considered 
had to be carefully selected for each crystal phase |13m 1139] . For the fluid phase we 
used 360 molecules. The number of molecules for each solid phase was chosen so as 
to fit at least twice the cutoff distance in each direction. Three different types of runs 
were performed: NpT, NVT Einstein crystal calculations and Gibbs Duhem integration. 
The equation of state (EOS) of the fluid was obtained from isotropic NpT runs, whereas 
anisotropic Monte Carlo simulations (Parrinello- Rahman like) [1241 11251 1126] were used 
for the sohd phases. In the NpT runs, about 40000 cycles were used to obtain averages, 
after an equilibration period of about 40000 cycles. However longer runs were used for 
the fluid phase at low temperatures. A cycle is defined as a trial move per particle plus 
a trial volume change. To evaluate the free energy of the solid, Einstein crystal (NVT) 
calculations were performed (with fixed center of mass). Also the length of the runs 
was of about 40000 cycles to obtain averages after an equilibration of 40000 cycles. In 
the Gibbs Duhem simulations a fourth order Runge-Kutta was used to integrate the 
Clapeyron equation. In total about 60000 cycles were used to pass from a coexistence 
point to the next one. When using Hamiltonian Gibbs Duhem integration 5-10 values 
of A where used to connect the initial to the final Hamiltonian. In the Gibbs Duhem 
simulations, the fluid, and cubic solid were studied with isotropic NpT runs whereas the 
solid phases were studied with anisotropic NpT runs. 

10.2. Free energy of liquid water 

The free energy of the liquid is calculated by integrating the free energy along a 
reversible path in which the water molecules are transformed into Lennard- Jones spheres 
by switching the charges off. The free energy of the reference Lennard- Jones fluid is 
reported in the work of Johnson et al. |220j . The energy (say, for the TIP4P model of 
water, being the treatment for SPC/E fully equivalent) of the system for a given point 
of the path. A, is given by: 

U{\) = XUTIP4P + (1 - X)Ulj (76) 

where A varies between (LJ) and 1 (water) along the integration path. Given that 
dA{X)/dX =< dU{X)/d\ 

>NVT^ the free energy difference between liquid water and the 
Lennard- Jones fluid is given by: 

Atip4p{N, V, T) - Alj{N, V, T) = {Utipap - ULj)My,T,x d^- (77) 
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where < Utipap — Ulj >n,v,t,x is an NVT simulation average computed for a given 
value of A. The integral is solved numerically (using Gauss-Legendre quadrature) by 
calculating the integrand at different values of A. In practice we perform the MC runs 
starting from A = 1 and going to A = 0. The final configuration of a run was used as 
the input configuration of the next run. The LJ fluid chosen as a reference state has 
the same LJ parameters {e/ks and a) as the water model. Therefore, the difference 
Utipap — Ulj is just the electrostatic energy. Once the Helmholtz free energy, A, is 
known, the chemical potential is obtained simply as: -jf^ = j, + J'^ j, . In this way 
we have computed the free energy of the liquid at 225 K and 443 K (see table [3]). In 
addition to the total free energy we report the free energy difference with respect to 
the reference LJ fluid (AA), the residual free energy for the LJ fluid and the ideal free 
energy {NkBT{hi{pA^) — 1)). In this work, for water, we shall assign the thermal de 
Broglie wave length to A = 1 A both for the liquid and for the solid phases. 

Table 3. Free energy of liquid water {Auquid) for the SPC/E and TIP4P models 
(with = qv = Qe = ^ )■ The residual and ideal contributions to the free energy 
of the reference LJ fluid are given. The residual term of the LJ fluid as obtained 
from the EOS of Johnson et al. [220] . The ideal term was obtained (in NkpT 
units) as ln(/jA^) — 1 where A = 1 A. The difference in free energy between the 
reference fluid and the water model AA is given. Simulations were performed in 
the NVT ensemble for the density d. 



Model 


T{K) 


p{bar) 


d{g/cm^) 


Aiiquid/iNkBT) 


Al^j/{NkBT) 


A^^j/iNksT) 


AA/{NkBT) 


SPC/E 


225 


564 


1.05 


-21.82 


2.500 


-4.350 


-19.97 


TIP4P 


225 


743 


1.05 


-19.48 


2.401 


-4.350 


-17.52 


SPC/E 


443 


4010 


1.05 


-9.53 


2.856 


-4.350 


-8.04 


TIP4P 


443 


4280 


1.05 


-8.58 


2.777 


-4.350 


-7.01 



Let us now present some consistency checks for the free energies. We shall only 
discuss it for the TIP4P model. We have computed the free energy at 225 K and 
443 K for the density d= 1.05 g/cm^. The free energy difference between both states 
is A443 K/{NkBT) - A225 K/{NkBT) = -8.58 + 19.48 = 10.90. Then we calculated 
the same difference by means of thermodynamic integration along an isochore (eq. [7]), 
obtaining again 10.90 NksT. Besides Jorgensen et al. have estimated the free energy for 
the TIP4P model at p=lbar and T=298 K to be G = -6.1 kcal/mol [22T]. Starting from 
the free energy at 225 K and d= 1.05 g/cw? and performing thermodynamic integration 
we obtained -6.09 kcal/mol, which is in excellent agreement with the value of Jorgensen 
et al. [22T] . 

Instead of using the LJ fluid, it is also possible to compute the free energy of the 
liquid taking the ideal g reference system. We obtain the free energy of TIP4P 

water at 240 K and 1.0174 g/cm^ by two different routes. In the first one the TIP4P 
is transformed into a LJ model. We obtained for the free energy of TIP4P in this state 
y4(240 K, 1.0174 g/cis?) /NkpT = -20.15. In the second route, a supercritical isotherm 
(T=900 K) is used from zero density to 1.0174 g/cm^ (obtaining using equation (12) 
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-4.699 for A/NksT of this intermediate state). Then we integrate along the isochore up 
to 240 K (with a free energy change computed by equation (7) of -15.434 NKbT). By 
adding these two numbers together we obtain A(240 K, 1.0174 g/cm^) /NkBT2 =-20.13 
from this second route, in very good agreement with that obtained by the first one. The 
computation of the free energy using the LJ fluid as a reference system is considerably 
shorter than using the ideal gas (besides this last route is especially difficult since the 
parameters of the Ewald sum should be chosen carefully along the supercritical isotherm 
integration). Values of the free energy of liquid water for other potential models have 
been reported recently [222] . 

10.3. Free energy of ice polymorphs 

The free energy of the different ice polymorphs is calculated using the Einstein crystal 
method (with fixed center of mass). The O is used as reference point and the field used 
is that described by equations (UHl), ([T9]) and (EOl). For disordered ices |223l ESI 1225] 
(Ih, Ic, III, IV, V, VI, VII, XII) the oxygens form a well defined lattice, but the water 
molecules can orient in different ways for a given oxygen lattice provided that the Bernal- 
Fowler ice rules |155j are satisfied. We have generated the disordered solid structures 
(with almost zero dipole moment) using the algorithm proposed by Buch et al. [226] 
(see Ref. |227j for other algorithm). For proton disordered ices we calculated the free 
energy for a particular proton disordered configuration. Due to the fact that there are 
many configurations compatible with a given oxygen lattice, there is a degenerational 
contribution to the free energy. The degenerational entropy of ice was estimated by 
Pauling in 1935 [156] as Sdeg = fcslnfi = ksN \n{3/2). Therefore, the disordered 
ice phases have an extra contribution to the free energy of —NkBThi{3/2). Ices III 
and V present partial proton ordering [228], and that decreases a little bit the Pauling 
estimate [2291 1230] . For ices III and V we shall use the sligthly lower value of the 
degenerational entropy estimated by MacDowell et al. [229j . For ices II, IX, VIII and 
the antiferroelectric analogous [231] of ice XI the protons are ordered, and there is 
no degenerational entropy contribution. Generating an initial configuration for proton 
ordered ices is relatively straightforward. 

The free energy calculations were performed in the NVT ensemble using the 
equilibrium shape of the simulation box obtained in anisotropic NpT runs. The location 
of the springs of the ideal Einstein crystal field were chosen to be close to the equilibrium 
positions/orientations of the molecules. The computed free energy should not depend on 
the particular choice of the positions and orientations of the ideal Einstein crystal field 
(provided that they are sufficiently close to the equilibrium position and orientations 
of the molecules in the absence of the external field). Several strategies are possible. 
For instance one could choose the position/orientations of the external field as those 
obtained from an energy minimisation at constant density (using the equilibrium box 
shape of the system). Another possibility is to use the experimental crystallographic 
positions (if available) of the atoms of the molecule and modify them slightly to satisfy 
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the bond lengths and angles of the model ( the bond lengths may be different in the 
model and in the real molecule, as for instance in SPC/E water). Also, experimental 
crystallographic positions (if available) could be used for the reference point (oxygen in 
the case of water) , and the orientations could be optimised from an energy minimization. 
We have used this last approach for ices. This approach has also been used recently 
by Wierzchowski and Monson for gas hydrates [68]. To obtain AA2 we used Gaussian 
integration with 12 points. AAi was evaluated from runs (200000 cycles) of the ideal 
Einstein crystal with fixed center of mass. The value of A^; (in ksT / {AY) was identical 
to the value selected for A^; ^ and A^; ^ (in /c^T units). The value of A^; was chosen so 
that the computed value of AAi differs by about 0.02NbT units from the lattice energy 
of the solid (defined as the intermolecular energy of the system when the molecules 
occupy the positions and orientations of the external Einstein field). 

Coexistence between phases at 150, 225 and 443 K were investigated. Therefore, 
the free energy calculations have been performed at those temperatures. In tables [Hand 
Owe report the free energy (Agoi/ (NkBT)) calculated for different ice phases for SPC/E 
and TIP4P respectively. The different contributions to the free energy [Aq, AAi, AA2) 
are given. The term Aq is the sum of A'^f^g_^^ plus AA^ plus a finite size correction 
(Frenkel-Ladd type, (2/A^)lnA^). For proton disordered ices the value of A^oi is the 
sum of Aq, AAi, AA2 and the Pauling degeneracy entropy —NksT ln{3/2). For proton 
ordered ices (XI,II,IX,VIII) the total value of A is just the sum of Aq, AAi, AA2. 
The lattice energy Uiatuce/NkBT is the energy of the solid when all water molecules 
remain fixed on the position/orientations of the Einstein crystal field (being the value 
of AAi/ (NksT) quite close to Uiattice/NksT). The free energies of the SPC/E model 
are lower than those of the TIP4P (this is consistent with their lower internal energies) . 

As a consistency check we determined via Einstein crystal calculations the free 
energy of ice VI at two different thermodynamic conditions for the SPC/E and TIP4P. 
Let us discuss the results for TIP4P. For TIP4P (ice VI) we obtained from free 
energy calculations Ai(225 K,1.480 g/cm3)=-17.67 NkeT and ^2(150 K,1.498 g/cm^)=- 

30.65 NksT (both states having a pressure of 25000bar). Starting from the value of 
Ai and performing thermodynamic integration we estimated ^2(150 K, 1.498 g/cm^)=- 

30.66 NksT in excellent agreement with the value obtained from Einstein crystal 
calculations. 

10.4- Determining the initial coexistence points 

Once the free energy of each phase is known, it is possible to find the points in the 
pressure-temperature plane at which two phases have the same chemical potential; i. 
e., the coexistence points. Given that most of the free energies and equations of state 
were obtained either at temperatures 150 or 225 K and at pressures 500 or 5000 bar, we 
focus the search of the coexistence points at these temperatures and pressures. 

Figure [5] shows the chemical potential as a function of the pressure at 150 K for 
ices Ih, II, and VI. For a given pressure, the phase of lowest chemical potential is the 
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Table 4. Free energy of the ice polymorphs for the SPC/E model (with qr = qv = 
Qe = 1). The free energy reported in the last column corresponds to the sum of all 
the terms {AQ+AA1+AA2) plus the degenerational free energy {—NkBTln{3/2)) 
for the case of orientationally disordered phases (the typical uncertainty of the 
resulting solid free energies is about 0.05NkBT). For ices III and V we did not 
use Pauling estimate for the degenerational free energy (— iVA;BTln(3/2)) since 
these ices present partial proton disorder, but the slightly lower value reported 
by MacDowell et a I. [229]. The number of molecules used for each solid phase is 
indicated in parenthesis just after the Roman numeral of the phase. A finite size 
correction (Frenkel Ladd type) has been included in ^o- The thermal de Broglie 
wave length A was set to A = 1 A. The residual internal energy of the ice U is 
reported, so that the entropy of the solid can be obtained easily from the relation 
S = {U — Asoi)/T. The orientational contribution to Aq (equation (29)) was 
computed from the approximate expresion given in Reference [48]. 
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most stable 


one. Ice Ih is 


the stable ph 


ase up 


to 3041 bar. 


At that pressure, ices 


Ih 



and II coexist. Beyond 3041 bar, ice II is the stable phase up to 6215 bar, where the 
chemical potentials of ices II and VI are equal. For higher pressures ice VI is the stable 
phase. By performing similar plots, coexistence points between different phases could be 
determined. In Table [6] these coexistence points are presented (for TIP4P and SPC/E). 
The relative stability between ices Ih and Ic (or between ices V and XII) could not 
be determined since the free energy difference between these solids was smaller than 
the typical uncertainty of our free energy calculations {0.05NkBT). As to the stability 
of ices Ih(Ic) with respect to ice XI (we used the antiferroelectric version of ice XI of 
Davidson and Morokuma |231j rather than the true ferroelectric version) , we found that 
the XI- Ih transition occurs at 84 K for SPC/E and 18 K for TIP4P (being the proton 
ordered ice XI the stable phase at low temperatures). 

The region of the phase diagram corresponding to 5000 bar (figure [6]) is the most 
problematic given that there are as many as seven phases competing; namely, ices II, 
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Table 5. Same as tabled for the TIP4P model. 



Ice 


p{bar) 


TiK) 


d{g/cm^) 


u 

NkuT 




An 


Ayl2 
NknT 


NkfiT 


-^sol 

Nkf{T 


Ih(288) 


1365 


150 
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Figure 5. Chemical potential as a function of pressure at 150 K for ices Ih, II and 
VI modeled with TIP4P. 

Ill, IV, V, IX and XII and the liquid. For SPC/E, ice II is clearly the most stable phase 
among the solid polymorphs. The liquid is again the stable phase at high temperatures 
(beyond 250 K). For the case of TIP4P, ices V and XII are the solid phases of lower 
chemical potential. The free energy difference between both polymorphs (0.008 NksT) 
is smaller that the error bar (0.05 NksT), so we could not determine which one is the 
most stable. Liquid water coexists either with ice XII or with ice V at 205 K at 5000 
bar. 
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Table 6. Coexistence points for TIP4P and SPC/E models from free energy 
calculations (and thermodynamic integration). 



Model 


Phases 


T(K) 


p (bar) 


SPC/E 


Ih-II 


150 


-444 


SPC/E 


II-VI 


150 


25270 


SPC/E 


liquid-II 


250 


5000 


SPC/E 


liquid- VI 


225 


20690 


SPC/E 


VI-VIII 


225 


57860 


SPC/E 


liquid-VIII 


225 


57500 


SPC/E 


liquid-VII 


443 


103520 


SPC/E 


liquid-Ill 
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500 


SPC/E 


liquid-Ic 


210 


500 


SPC/E 


liquid-XI 


187 


500 


TIP4P 


liquid-Ih 


228.8 


500 


TIP4P 


Ih-II 


150 


3041 


TIP4P 


II-VI 


150 


6215 


TIP4P 


II-V 


152.6 


5000 


TIP4P 


II-III 


180.3 


5000 


TIP4P 


liquid-Ill 
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5000 


TIP4P 


liquid-V 


204.1 


5000 


TIP4P 


liquid- VI 


225 


8940 


TIP4P 


VI-VIII 


225 


57290 


TIP4P 


liquid-VII 


443 


91940 


TIP4P 


liquid-Ic 


228.8 


500 


TIP4P 


liquid-XI 


192 


500 


TIP4P 


liquid-XII 


205.0 


5000 



10.5. The phase diagram of water 

Once an initial coexistence point has been determined, by using Gibbs Duhem 
integration (either dp/dT or dT/dp ) it is then possible to draw the complete phase 
diagram. In certain cases the point where two coexistence lines met (triple point) 
was used as origin of the third coexistence line emerging from the triple point. The 
complete phase diagram of SPC/E and TIP4P is presented in figure [71 As can be 
seen SPC/E fails in reproducing the phase diagram of water (notice for instance that 
ice Ih is stable for this model only at negative pressures), whereas TIP4P provides 
a qualitatively correct description of the phase diagram (except for the high pressure 
region of the phase diagram ). The reason of the failure of SPC/E and success of TIP4P 
has been identified recently. The low quadrupole moment of SPC/E and the high value 
of the ratio dipole/quadrupole of this model is the cause of the failure |112[lll4j . In fact 
TIP4P provides a quadrupole moment and a ratio dipole/quadrupole in much better 




T(K) T(K) 

Figure 6. Gibbs free energy versus temperature at 5000 bar for different phases of 
SPC/E (left) and TIP4P (right). 
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Figure 7. Phase diagram of water as obtained from experiment (center) and from 
computer simulations for the TIP4P model (left) and for the SPC/E model (right). 
The filled circles on the right panel indicate the stability limit of the solid phases in 
NpT simulations (without interfaces). Notice the shift of 100 MPa in the right panel. 

agreement with experiment. The effect of a quadrupole moment on the vapor-hquid 
equihbria of molecular models is well known |232l I233[ 1234] . However it seems that 
the role of the quadrupole on water properties has been overlooked in spite of some 
warnings about its importance [2351 1236i I237i I218i I238i I239j . Figure [7] illustrates how 
the evaluation of the phase diagram of water by computer simulation is indeed possible. 

10.6. Hamiltonian Gibbs Duhem simulations for water 



The liquid-Ih solid coexistence temperatures at p = 1 bar for TIP4P and SPC/E 
have been estimated from free energy calculations to be T = 232 ± 5 K and T = 
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215 ± 5 K, respectively. These numbers are in relatively good agreement with estimates 
from other authors for TIP4P [221 1231 [W] and for SPC/E [HSl UHl flOT] . 
An interesting question is whether these two values (for TIP4P and SPC/E) are 
mutually consistent. Starting from the SPC/E model and performing constant pressure 
Hamiltonian Gibbs Duhem simulations (integrating the generalised Clapeyron equation 
as described previously) one should recover the melting temperature of the TIP4P. In 
fact starting from the SPC/E ice Ih melting point we obtain T = 232.3 K for TIP4P (see 
figJHl ) in very good agreement with the result obtained through free energy calculations. 




• # From SPC/E to TIP5P 

♦ ♦From TIP4Pto TIP5P 



Figure 8. Hamiltonian Gibbs Duhem integration results. Left. Melting temperature 
of ice Ih as a function of the parameter A connecting the models SPC/E (A = 0) 
and TIP4P (A = 1). The points were obtained by using Hamiltonian Gibbs Duhem 
integration. The dashed line is a guide the eye. The horizontal lines correspond to 
the melting temperatures of SPC/E (dotted line) and TIP4P (dashed-dotted line) as 
obtained from free energy calculations. Right: Melting temperature of ice Ih for the 
TIP5P model (A = 1) obtained from Hamiltonian Gibbs Duhem integration starting 
from SPC/E or TIP4P models. When connecting two water models by Hamiltonian 
Gibbs Duhem integration, the position of the oxygen atom and of the HOH bisector 
was the same in both models. 



Once the melting point of ice Ih for TIP4P and SPC/E seems to be firmly 
established one could use these values to estimate (by using Hamiltonian Gibbs Duhem 
simulations) the melting point of another water model, as for instance TIP5P |241j . 
Obviously, the properties of the final model should be independent of the reference 
model. When the starting model is SPC/E we obtain T = 275 K for TIP5P whereas 
the calculated result using the TIP4P model as a reference is T = 273 K. The agreement 
between both estimates is satisfactory taking into account that the error of the Gibbs- 
Duhem integration is about 3 K. This is is further illustrated in the right panel of 
figure [8] which shows the results of the integration. By using Hamiltonian Gibbs Duhem 
integration the melting point of ice Ih for other models of water was determined. They 
are presented in Table [71 Notice that most of the water models tend to give low melting 
points. 

For models with three charges (SPG, SPC/E, TIP3P, TIP4P, TIP4P/Ew, 
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Table 7. Melting points obtained from free energy calculations (TIP4P and SPC/E), 
Hamiltonian Gibbs Duhem integration (rest of the models) [215| I211|, 12421 1214j , and 
from direct fluid-solid coexistence. The last column is the value of Ts (see text) 
obtained from simulations of ice Ih with a free interface. TIP4P-Ew is a water model 
proposed by Horn et al. |243j and NvdE is a six sites model proposed by Nada and 
van der Eerden [244j . 



Model Free energy Direct coexistence Free surface 



TIP4P/Ice 


272(6) 


268(2) 


271(1) 


TIP4P/2005 


252(6) 


249(2) 


249(3) 


TIP4P-EW 


245.5(6) 


242(2) 


243(2) 


TIP4P 


232(4) 


229(2) 


230(2) 


TIP5P 


274(6) 


271(2) 




SPC/E 


215(4) 


213(2) 


217(2) 


NvdE 


290(3) 


288(3) 






Figure 9. Left: melting temperature of ice Ih at 1 bar plotted as a function of the 
quadrupole moment Qt (taken from )I12j ). Qj- is defined as the value of {Qxx~Qyy) 1'^ 
of the quadrupolar tensor, where the x axis joins the two H atoms and the y axis is 
perpendicular to the molecular plane. Right: density of water at room pressure as a 
function of temperature as obtained from experiment (filled circles) and from computer 
simulations of several water models (lines). The open triangles indicate the melting 
point of ice Ih for each model. 



TIP4P/Ice, TIP4P/2005) a correlation between the melting point and the quadrupole 
moment has been found. This is illustrated in figure [9l It is seen that models with rather 
low quadrupole moment (TIP3P, SPC, SPC/E) provide rather low melting points. The 
melting temperature of TIP4P is closer to the experimental value. Motivated by this we 
have proposed a new modified TIP4P model, with a higher quadrupole moment, able to 
reproduce the experimental melting point of water. We have denoted this new model as 
TIP4P/Ice |214j . A second finding was that for three charge models, the temperature at 
which the maximum in density occurs at room pressure (TMD) is about 20-25 K above 
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the melting temperature |242l 1245] . Experimentally for water, the maximum in density 
occurs 4 degrees above the melting point. Therefore it is impossible with three charge 
models to reproduce simultaneously the melting point and the TMD. It is likely that the 
inclusion of quantum effects and/or polarizabihty [MSI IMS [2501 ISMl [252] may 

be needed to reproduce these two properties simultaneously. For this reason we have 
proposed the TIP4P/2005, which reproduces the TMD of real water better than any 
other water model proposed so far (see figureQb). An interesting question is to analyze 
whether these new models still predict correctly the phase diagram of water. By using 
Hamiltonian Gibbs Duhem integration it is possible to estimate the phase diagram of 
a certain water model, starting from the phase diagram of another reference model. 
Thus by using TIP4P as reference, we have estimated the complete phase diagrams 
for TIP4P/Ice and TIP4P/2005. The obtained phase diagrams are presented in figure 
[TOl As can be seen these models predict quite well the fluid-solid equilibria of water 
improving the predictions of TIP4P. The TIP4P/2005 yields also an excellent prediction 
of the vapor-liquid equilibria. 




Figure 10. Phase diagram for the TIP4P family, a) Left paneh fluid-sohd equihbria 
(the hnes show the predictions for several TIP4P-like models and the symbols represent 
the experimental data), b) Right panel: vapor- liquid equilibria (the symbols show the 
predictions for TIP4P-like models and the lines represent the experimental data). 



10.7. Direct coexistence simulations 

To estimate the melting point of ice Ih for several water models by direct coexistence, 
NpT simulations were performed with 870 molecules and the MD program Gromacs 
|253[ I254j . In the initial configuration half of them formed ice, and the other half were 
in the liquid state. Both phases were in contact so that these are direct coexistence 
simulations. 

The evolution of the energy for the TIP4P/2005 model with time is presented in 
figJTTl As can be seen the energy increases with time for T = 252, 254, 256 K reaching a 
plateau (the plateau indicates the complete melting of the ice) . The energy decreases for 
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Figure 11. a) Left panel. Evolution of the total energy (per mole of molecules) 
in NpT MD simulations of a box containing ice and liquid water at 1 bar for the 
TIP4P/2005 model, b) Right panel. Total energy as a function of time obtained at 
several temperatures by performing MD simulations of TIP4P/2005 for a block of ice 
Ih with a free surface. 



T = 242 K reaching a plateau (the plateau indicates the complete freezing of the water) . 
Snapshots of these final configurations can be found in Ref. |189j . At a temperature 
of 249 K the energy does not change with time, and the interface is stable after 10ns. 
Therefore this is the estimate of the melting point by direct coexistence for TIP4P/2005. 
Similar runs were performed for other water models. In Table the melting points 
of different water models as estimated from direct coexistence are compared to those 
obtained from free energy calculations. The agreement between both techniques is quite 
good. Direct interface simulations have been used by several authors to estimate melting 
points or ice growth rates for different water models [921 1255[ I194[ I256[ I187[ I192j . 

10.8. Melting point as estimated from simulations of the free surface 

In figure [TT] the evolution of the total energy of ice Ih (having a free surface) with 
time is presented |209j . At high temperatures, the total energy of the system increases 
continuously and then reaches a plateau (that corresponds to the complete melting of 
the solid). The behaviour at low temperatures is different. At the beginning (first 1- 
2ns), there is an increase of the energy but after that the energy remains approximately 
constant, apart from the thermal fiuctuations. The analysis of the configurations of the 
TIP4P/2005 at T = 245 K, shows that the increase of energy during the first Ins is due 
to the formation of a thin liquid layer at the surface of ice, which may indicate the onset 
of surface melting, mentioned already in the Introduction, and first proposed by Faraday 
|257] . The formation of a quasi liquid layer on the surface of ice below has been found 
both in experiment (see [258[ I259[ 12601 126H I262j and references therein) and in computer 
computer simulation for several potential models of water [2631 1264[ I265i I266[ [92] and 
it has been explained by several theoretical treatments [2581 I260j . By repeating the 
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simulation at several temperatures it is possible to determine the lowest temperature 
at which the block of ice melts T+, and the highest temperature at which it does not 
melt T_. By taking the average of these two temperatures we obtain what we call 
Tg = (T_|_ + T_)/2. Tg provides an estimate of the melting point. The values of Tg 
obtained for water models are presented in Table [71 As can be seen, Tg is identical 
to Tm, within the error bar. Thus for ice Ih, the presence of a free surface suppresses 
superheating and ice melts at the equilibrium melting temperature (although runs of 
about 10ns or longer may be needed). In figure [I2] the final configuration (after a 8 ns 




Figure 12. Instantaneous configuration of the TIP4P/Ice system at T = 268 K at 
the end of a 8ns run. Although the temperature is well below the melting point of the 
model, a quasi-liquid layer is clearly present in the ice-vacuum interface. 

run) obtained for the TIP4P/Ice at a temperature well below the melting point of the 
model (T = 264 K). As can be seen, a quasi liquid layer is already present in the system. 

10.9. Properties at K 



Table 8. Residual internal energy (in Kcal/mol) of several ice polymorphs atT = K 
and p = for popular water models. The results for the most stable phase of each 
model are presented in bold. 



Ice 


TIP4P/Ice 


TIP4P/2005 


SPC/E 


TIP5P 


Ih 


-16.465 


-15.059 


-14.691 


-14.128 


II 


-16.268 


-14.847 


-14.854 


-14.162 


III 


-16.140 


-14.741 


-14.348 


-13.320 


V 


-16.049 


-14.644 


-14.169 


-13.101 


VI 


-15.917 


-14.513 


-13.946 


-12.859 



In Table[8]the residual internal energies at zero T and p are given for the TIP4P/Ice, 
TIP4P/2005, SPC/E and TIP5P models. For TIP4P/Ice and TIP4P/2005 ice Ih is the 
structure with the lowest energy (probably ice XI which is a proton ordered version 
of ice Ih would have an slightly lower energy but it was not considered for this study). 
However for SPC/E and TIP5P the lowest internal energy corresponds to ice II. Thus, for 
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TIP4P/2005 and TIP4P/Ice ice Ih is the stable phase at zero pressure and temperature 
whereas for SPC/E and TIP5P the stable phase is ice II. From the properties at K 
it is simple to locate the Ih-II transition pressure at K (see |267] ). In figure fT3l the 
predicted pressures using the calculations at K are presented (circles). The lines 
represent the results obtained from free energy calculations (used to obtain an initial 
coexistence point at 150 K) and Gibbs Duhem integration. It may be seen that both 
sets of calculations agree quite well so both sets of results are mutually consistent (i.e., 
the estimated coexistence pressure at K is the same). It is clear that ice II is more 
stable than ice Ih at zero temperature and pressure for the SPC/E and TIP5P models. 
This example illustrates how K properties can be used to test for self consistency in 
phase diagram predictions. They are also quite useful to test the performance of water 
models [26HI [HHl EQI ESSl l267j . 
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Figure 13. Coexistence lines between ices Ih and II as obtained from Gibbs-Duhcm 
simulations for TIP4P/2005, TIP4P/Ice, SPC/E and TIP5P models (solid lines). The 
symbols represent the coexistence pressures as obtained from the properties of the 
systems at zero temperature. For each water model, ice Ih is the stable phase below 
the coexistence line (low pressures) and ice II is the stable phase above the coexistence 
line (high pressures). 




11. Phase diagram for a primitive model of electrolyte 

Let us now present another example of a phase diagram calculation for a completely 
different model, the restricted primitive model (RPM). The restricted primitive model 
(RPM) is one the simplest model of electrolytes. It this model the cations are represented 
by hard spheres of diameter a having a charge +g and anions represented by hard 
spheres of diameter a having a charge —q. The model is quite simple and for this reason 
it can also be studied theoretically [2701 1271] . The system has vapor-liquid equilibrium 
|272[ 12731 1274j (in spite of the absence of dispersive forces). Several solid structures can 
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be considered |275] . the simplest being the CsCl hke structure (a bcc type of structure 
with the anions occupying the vertices of a cube and the cations occupying the center of 
the unit cell). Another possible structure is the fee disordered structure. In this structure 
the ions occupy an fee lattice, but with positional disorder (cations and anions occupy 
the lattice points in a disordered way). At low temperatures it is possible to conceive 
a solid structure with an fee arrangements of the ions, but with positional order. The 
symmetry of the phase is tetragonal. This phase will be labelled as the tetragonal 
phase [19]. It is presented in figure [TH Free energy calculations (Einstein crystal) were 
performed to determine the free energy of the CsCl and tetragonal structures. Due to 
the presence of partial disorder the Einstein crystal method can not be applied directly 
(to a snapshot) to get the free energy of the fee disordered structure (although one may 
suspect that an strategy similar in spirit to that proposed for the plastic crystal phases 
can be also successful here if the external field is able to lead the system from disordered 
configurations to an ordered solid without crossing first order transitions). The RPM 
system becomes a hard sphere at infinitely high T for which the free energy is known 
(a trivial mixing contribution should be added). For this reason the free energy of the 
fluid and of the fee disordered solid can be obtained by thermodynamic integration. 
Exchange moves (where a cation and an anion exchange their positions) were used to 
sample correctly the disorder, both in the fluid and in the fee disordered structure. 




Figure 14. Left: Tetragonal ordered structure of the RPM model. The atoms 
form an fee structure, but the ions are ordered. Right: Phase diagram of the RPM 
model showing the equilibria between vapor and liquid (inverted triangles), fluid and 
CsCl like structure (filled circles), CsCl like structure and tetragonal (fee ordered) 
phase (rhombs), fluid and tetragonal structure (open circles), fluid and fee disordered 
structure (triangles), and ordered-disordered fee phases (squares). 

After computing the free energies, some initial coexistence points were determined, 
and then by using Gibbs Duhem integration the complete phase diagram was computed. 
The resulting phase diagram |276i I277i H9| [50] is presented in figure [HI At 
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high temperatures, freezing leads to the substitutionally disordered close packed 
structure. By decreasing the temperature this solid structure undergoes an order- 
disorder transition transforming into the tetragonal solid. At low temperatures freezing 
leads to the caesium chloride structure (CsCl) which undergoes a phase transition to 
the tetragonal structure at high pressures. The tetragonal solid is the stable sohd phase 
at low temperatures and high densities. In a narrow range of temperatures coexistence 
between the fluid and the tetragonal solid is observed. Three triple points are found for 
the model considered: the usual vapor-liquid-CsCl, the fluid-CsCl-tetragonal and the 
fluid-fcc disordered-tetragonal triple point (notice on figure [T3] the narrow range of the 
fluid- tetragonal solid coexistence line). 

Although initially conceived to describe electrolytes and ionic salts, the RPM 
has been found to be quite useful to describe certain colloidal mixtures, consisting 
in equimolar mixtures of colloidal particles of equal size, but with charges of different 
sign. Thus a colloidal version of the RPM exists |278] . For these colloidal mixtures, 
three different solid phases have been found experimentally, the fee disordered solid, the 
CsCl solid, and a fee ordered structure |279l 1280] . The fee ordered structure was found 
to be of CuAu type, rather than the tetragonal structure of figure [TH This difference 
with the RPM phase diagram seems to be due to the fact that in colloidal mixtures 
the interaction between charged particles is of Yukawa type rather than being purely 
Coulombic. This is so because the interaction between charged particles is screened by 
the presence of an ionic atmosphere due to the counter ions of the colloids. In fact, 
when this screening is incorporated in the potential with a Yukawa type model, the 
CuAu structure was found to be stable in a thermodynamic region between the CsCl 
and the tetragonal structure [56l I281j . 

In summary the RPM has proved to be quite useful in the description of mixtures 
of charged colloidal particles. It would be of interest to determine the phase diagram 
for a model where the two spheres of the model present different size. This is usually 
called the primitive model. In fact the primitive model (where cations and anions have 
different size) may indeed be a more general model than RPM, since in real ionic fluids 
the size of the ions is usually different. A colloidal realization of the PM model has been 
obtained experimentally |279j . Finally, there is an increasing interest in determining the 
properties of ionic liquids. It would be of interest to determine the factors affecting the 
melting point of ionic liquids |282j which are regarded as new solvents. |283j . Work in 
this direction has appeared recently |284] . 

12. Phase diagram of a simplified model of globular proteins 

A flnal example of the application of the techniques described here is provided by the 
calculation of the phase diagram of a simple model of globular proteins. During the 
last few years there has been an increasing number of studies of the phase behaviour of 
globular proteins using very simplifled models. The flrst studies have been performed 
using short ranged isotropic potentials and, even with these very simple models, it was 
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already possible to reproduce some of the features of the phase diagram of proteins, e.g., 
the existence of a metastable critical point |285j . However, proteins are known to form 
very low density crystals, with densities below those of the close-packed crystals typically 
formed by isotropic potentials, which is indicative of highly directional interactions [286j . 
Further evidence of the importance of anisotropy in the interactions among proteins has 
been recently obtained in a theoretical study of the fiuid-fiuid equilibria. This study has 
shown that a quantitative agreement with experiments is obtained by the introduction 
of anisotropy [287] . as opposition to isotropic models that only provided a qualitative 
description. Moreover, theoretical studies of the phase behaviour of anisotropic models 
are also acquiring much interest due to the fact that several experimental groups have 
recently been able to produce colloids that are anisotropic either in shape or in their 
interactions |288l 12891 12901 1291j . So far there has been already a few theoretical studies 
of the phase behaviour of simple anisotropic models [IHl [2921 [2931 |295l [581 [2961 [297] , 
although most of them were concerned with the fluid-fluid equilibria rather than with 
the solid-fluid and solid-solid equilibria. 

We have used a very simplified model which consists of a repulsive core (the 
Lennard- Jones repulsive core), plus an attractive tail modulated by Gaussian functions 
located at some given positions or patches, which will be specified by some vectors. The 
total energy between two interacting particles will be given by [298[ I299[ [57] : 



where Vlj is the Lennard- Jones potential, (Jpatchy is the standard deviation of the 
Gaussian, 6k,ij iOiji) is the angle formed between patches k {I) on atom i (j) and the 
interparticle vector Yij ( r^j), and kmin i}min) is the patch that minimises the magnitude 
of this angle. The interaction is a maximum when both patches are pointing at each 
other along the interparticle vector rjj and it will decrease as the particles deviate 
further from this equilibrium orientation. We have chosen to study a model with 
6 patches distributed in an octahedral symmetry. A relatively narrow width of the 
patches (crpatc/ty=0.3 radians), for which it is expected that the low density simple cubic 
(sc) crystal becomes stable. 

Besides the sc crystal in which each of a particle's patch is pointing at each one of 
its six nearest neighbours (see figure [T5l) . there are other three solid phases that might 
be formed with this model and at this patch width, (Jpatchy = 0.3 radians. The first 
structure is a body centred cubic (bcc) solid, in which each patch is aligned with the 
second neighbours. This structure can be also seen as two interpenetrated sc lattices 
that almost do not interact between each other (similar to the behavior of high density 
ice polymorphs). Therefore, a higher density crystal is obtained with a low penalty in 
the energy. At high pressures it is expected that a closed packed face centred cubic 
(fcco) structure will also appear. In this case the patches will be also pointing to the 




Tij < (TlJ 
Tij > (TlJ 



(78) 




(79) 
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second neighbours. This structure exhibits a much higher energy than both the sc and 
the bcc sohds. Finally, at high temperatures, it is expected that a plastic phase will 
also appear (feed), i.e., a solid where the center of mass of each particle is located at 
the lattice positions of a fee structure, but that is orientationally disordered. 
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Figure 15. Left: Phase diagram of the octahedral six-patches model. Labels show 
the region of stability of each phase. The dashed line in the fluid-sc coexistence 
region signals the expected fluid- fluid binodal (see Ref. [S^ for more details). Right: 

Structure of the low density sc solid for the octahedral anisotropic model. 



Following the procedure described before, the coexistence point between two phases 
was computed by imposing the conditions of equal temperature, pressure and chemical 
potential. The free energy of the fluid was computed by thermodynamic integration. 
In the case of the fluid, the equation of state was integrated up to very low densities, 
where the fluid can be considered to behave as an ideal gas. The free energy of the solid 
was computed by Einstein crystal calculations. Once that a coexistence point is known 
the whole coexistence line have been traced using the Gibbs-Duhem method. 

The resulting phase diagram is plotted in figure [151 At high temperatures the fluid 
freezes into the orientationally disordered plastic crystal phase (feed), at intermediate 
temperatures into the bcc solid and at low temperatures into the sc crystal. The sc 
structure is destabilised at high pressure by the bcc solid and, at even higher pressures, 
the ordered fee solid becomes stable (fcco). This fcco ordered solid transforms into 
a plastic crystal feed as the temperature increases. There are three triple points in 
the phase diagram: the fluid-sc-bcc, fluid-fccd-bcc and bcc-fccd-fcco. Finally, it is also 
worth noting that in the neighbourhood of the fluid-sc-bcc and fluid-fccd-bcc triple 
points, reentrant melting occurs. Coexistence points from free energy calculations were 
found to be in agreement with those found from direct coexistence calculations. 

In summary, even for a relatively simple model potential, we have obtained a 
complex phase diagram with many solid phases and unusual behaviour such as reentrant 
melting. But the most interesting finding is that even with a very simple model as the 
one described here, it is possible to reproduce two important features of the phase 
diagram of globular proteins, namely, the existence of a metastable critical point and 
the stabilisation of a low density crystal. Similar behavior has been found for a primitive 
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model of water |146l HSl I3UU] . Studies of nucleation of these models can be very relevant 
to understand the crystallization of globular proteins. 

13. Conclusions 

In this work we have reviewed the methodology of free energy calculations and 
given examples of recent work on the determination of fluid-solid and solid-solid 
equilibria by computer simulations. Free energy calculations are used to compute 
initial coexistence points between phases, and then Gibbs Duhem integration is used to 
compute coexistence lines. Other procedures to estimate fluid-solid equilibria as direct 
coexistence and simulations of the free surface of the solid have been discussed. The 
Einstein crystal methodology and the Einstein molecule approach have been presented 
in a rather comprehensive way. Both methodologies yield identical values of the free 
energy. It is shown that the free energy presents a strong N dependence, and that finite 
size corrections are needed to estimate properties of solids in the thermodynamic limit. 
The issue of the symmetry of the orientational field in Einstein crystal calculations for 
molecular fluids has been discussed. We do hope that the extensive discussion of all 
these aspects helps other researchers in the area to perform free energy calculations and 
phase diagram determination. In fact there are at least six areas where one can predict 
intense activity in the future. The flrst one is the determination of the phase diagram 
of molecular fluids. In this work, the procedure to obtain free energies for water is 
presented. Besides, free energies and coexistence points for SPC/E and TIP4P models 
of water are presented by the first time. These results lead to the determination of the 
full phase diagram for water, performed recently by our group [IIOJ. The example of 
water illustrates clearly that phase diagram calculations for molecular fiuids is indeed 
feasible nowadays and that it can help to improve current models. In fact these free 
energy calculations led to the proposal of an improved version of TIP4P denoted as 
TIP4P/2005. This model is able to describe correctly the phase diagram of water, 
the maximum in density of water, the density of the ice polymorphs including the 
methane hydrate |301l 13021 1303] , the vapor-liquid equilibria [304] , the surface tension 
[305] ■ the diffusion coefficient, and the structure of water |215l I306j over a wide range 
of temperatures and pressures. The determination of the phase diagram of TIP4P was 
a crucial step in the development of TIP4P/2005. We do not see any difficulty in 
performing similar studies to improve potential models for other molecular fluids. The 
work on water, shows that melting points obtained from free energy calculations, direct 
coexistence simulations and free surface simulations, are almost identical (taking into 
account the statistical uncertainties and the slightly different implementations of the 
potential). The second area where phase diagram calculations can be useful is in the 
study of ionic systems. Here we reviewed the phase diagram of the RPM model, but it 
is clear that it would be of great interest the determination of the phase diagram of PM 
models and of other models of salts (including probably ionic liquids which are becoming 
increasingly important from a technological point of view). Some ionic systems can also 
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be used to describe charged colloids. 

The third area is the area of crystallization of proteins. It is now clear that models 
with short range anisotropic forces can be regarded as primitive models of proteins. 
In such models the liquid-liquid separation is metastable with respect to freezing, and 
the competition between phase separation and crystallization is of great interest to 
understand the presence or absence of crystallization in proteins. The fourth area is 
the study of freezing under confinement due to the interest in understanding fluid- 
solid equilibria on the nanometer scale |307[ I308[ 1309] . The fifth area is the study of 
the solubility of salts(or solids in general) in water (or solvents in general) where the 
knowledge of the chemical potential of the solute in the solid phase is required. Very 
little effort has been devoted to this problem |143i I147j . Finally studies on nucleation 
[120] should be benefit from the knowledge of the equilibrium melting temperatures. In 
summary, the study of fluid-solid and solid-solid equilibria of molecular and complex 
systems by computer simulation is now feasible, and the procedures to do it seem well 
established. The study of fluid-solid and solid-sohd equilibria by computer simulation 
can play a central role in developing potential models for condensed phases and for 
providing molecular understanding of a number of phenomena involving solid and liquid 
phases. The enormous sensitivy of phase diagrams to interaction potentials allows to 
test the performance of the different potentials available for a certain substance, and 
offers a unique opportunity for their improvement. 



14. Appendices 

14- 1- Appendix A. Partition function of the Einstein crystal with fixed center of mass 



The translational contribution to the partition function of an Einstein crystal with fixed 
center of mass is: 
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The integral over the space of momenta is not relevant to compute the free energy and, 
therefore, we will leave aside this contribution, that we will include simply as a factor 

pCM. 
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We will focus on the integral over the configurational space: 
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This integral can be expressed in a simpler way by defining a change of variable, 
— i"io = r[. The Jacobian of this change of variable is 1, and the configurational 
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The Dirac delta function can be written [310] as: 
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the configurational integral can be written as: 
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Each term in the summation can be rewritten in a more convenient form: 
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Notice that in equations (84-86) there is an implicit scalar product between the vector 
k and the accompanying vector. The integral can then be expressed as: 
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It is now convenient to do another change of variable: 
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To compute the Jacobian associated to this transformation, we will consider the simple 
case of a system in one-dimension and with two particles. The Jacobian J is given by 
the following determinant: 
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After applying the change of variable, the integral can be expressed as: 
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This integral can be split in two Gaussian integrals: 
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whose solution is: 

m/2 / .n f. \ 3/2 



TT 



- (2^)3 Ua J lEfir/i' ^^^^ 



Doing a bit of algebra it can be shown that: 

3/2 / ^ \ N3/2 / AT \ "3/2 

TT J 



/ iV 
\i=l 



z£i = F^ hM (E^"?! (95) 



or, more simply: 

/ ^ \ 3(iV-l)/2 / AT X -3/2 

Summarising, we have obtained that the partition function of an Einstein crystal with 
fixed center of mass is given by: 



Q%rn,=p''''(^] (E/^n • (97) 



3(Af-l)/2 / AT s -3/2 



When all molecules are identical the reduced mass fii is simply Therefore the 

previous equation can be simplified to: 

/ ^ \3(7V-l)/2 

Qif^, = P'''' j (Nf' . (98) 

which is the final expression for the free energy of an ideal Einstein crystal with fixed 
center of mass. An explicit expression for P*^*^ is not needed to get the free energy of 
the solid since it cancels out with a similar term in equation fH6|l . However, it is not 
difficult to obtain P^*^ by realizing that equation flSTl) is formally identical to equation 
(!82|) (with Hi = 1 and A^; = l/(2mj) and ommiting the prefactor h^^^~^^). A derivation 
similar to that used to get equation (!96|) from equation (182!) leads to: 

P''"' = ^N'"' (99) 
to be compared with 

^ = (100) 

so that the ratio p^^^ / p adopts the value A^A^~^/^. If equation (p9|) is replaced into 
equation (l98il one obtains: 

3{Ar-l)/2 

A3(^-i) \i3A~ 



«. = T4zTTh^l (101) 



which is dimensionless as it should be. 
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14-2. Appendix B: Computing within a Monte Carlo program 

The condition of fixing the center of mass of the reference points is quite useful since 
it ehminates any divergence of the integrand of equation ( l39i) . A displacement c? of a 
given molecule must be accompanied by a displacement {—d/{N)) of all the molecules 
of the system (assuming all particles are identical), so that the center of mass remains 
in its original position. In practice, this is not very convenient. It is more convenient to 
perform a simulation without the restriction over the center of mass but keeping track 
of the position of the center of mass [139J. Let us denote as the initial position of the 
reference point of molecule i in the perfect lattice, and AUcm is the difference between 
the present position of the center of mass and its initial position. Let us denote by rf 
the actual position in the simulation (without the restriction in the center of mass) of 
the reference atom of molecule i. Let us define Ar, = rf — r^o — AUcm- Let us compute 
the energy change when in a trial move, the random displacement of molecule is Aj. 
The energy of the old (prefix old) and new (prefix new) configurations is given by: 

U^T'/^E =(Arf^)2 + E(Arf)^ 

U^^nrV^E = (rf + A. - r.. - ARf,, - ^) ^ + E (r^'" - r,o - AR^t, - 

= [^rf + A. - + E (Arf - (102) 

We have assumed that all particles have the same mass, so that a displacement Aj of 
one of them leads to a displacement Aj/A^ of the center of mass. We have not included 
the orientational contribution since it cancels out when computing the energy change 
(i.e., the orientational energy is not affected by the translation of molecule i). Therefore, 
the potential energy change will be given by: 




2Arf ^ 



(103) 



This equation can be simplified since, as the center of mass is constrained, it holds that 
Y.f=i Arj = (in the right hand side, the second and last terms cancel out): 

AUEtl/AE = 2Arf . A. + ^ + ^ (^) ^ ' 

It is easy to show that this expression can be also written as: 



AUEtl/^E = 2Ar,f . A, + A,^ — — (105) 



/N -V 

In this way, it is possible to perform a MC simulation without keeping the center of 
mass fixed, but including this constraint through the Monte Carlo acceptance rule. 
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14-3. Appendix C. The Frenkel-Ladd expression 

In 1984 Frenkel et Ladd derived an expression for the free energy of a solid. It is not the 
same as that given by Poison et al. |142] and presented also in this paper. The reason 
of this difference is the following: 

• 1. The expression used for AA3 by Frenkel and Ladd was (l/A^)lnA^ (in NksT 
units) lower than the correct one. 

• 2. The expression used for A'^f^_^^ by Frenkel and Ladd was (3/A^) In (in NksT 
units) higher than the correct one. 

Taking into account both contributions it turns out that the Frenkel Ladd expression 
gives an energy (in NksT units) {2/N) InN higher than the correct one. Let us describe 
briefly the source of these two discrepancies. For AA^ Frenkel and Ladd used : 

A^f^ = Asoi - Agf = A;BT(ln(P^*VP) - ln{V)) (106) 

instead of : 

AA, = Asoi - AZY = kBT{HP^''/P) - \n{y/N)) (107) 

which is the expression to be used when all permutations, A^!, are included in the 
reference ideal Einstein crystal. The second discrepancy is due to the fact that the 
constraint of fixing the center of mass was implemented by Frenkel and Ladd as : 

N 

^(r, - r,,) = (108) 

i=l 

instead of : 

N 

^/i,(r,-rio) = (109) 

i=l 

with /ij = 1/A^. One can simply say that to fix the center of mass Frenkel and Ladd 
used yUj = 1 instead of fii = 1/N. It is simple to analyze the mathematical consequences 
of that. It is just enough to look Appendix A, and to see what happens in the final 
expression (equation (!97|) ) when /ij = 1 is used instead of fii = 1 / N . Then one obtains: 

/ ^ \3(iV-l)/2 

QErn, = P''"\-^J (iV)-^/^ (110) 

By comparing equation (11 101) (Frenkel-Ladd) with equation (l98l) it is simple to see how 
the Frenkel Ladd expression for Q^fnt gives a contribution (in NksT units) (3/A^) InA^ 
higher than the correct one. 

In summary, when all factors are considered, the Frenkel-Ladd expression gives an 
energy (in NksT units) {2/N) InN higher than the correct one. 
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